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ABSTRACT 


Advanced  techniques  have  been  employed  for  decades  in  the  location  of  emitted 
communications  signals.  The  use  of  cyclostationary  maximum-likelihood  estimators 
(MLE)  applied  to  time  difference  of  arrival  (TDOA)  and  frequency  difference  of  arrival 
(FDOA)  techniques  for  target  emitter  location  is  a  newly  developed  approach.  This 
dissertation  provides  an  introduction  to  traditional  TDOA  and  FDOA  methods  and  the 
underlying  cyclostationary  theory  of  pseudo-random  signals.  A  survey  of  current 
cyclostationary  TDOA  and  FDOA  methods,  the  description  of  new  maximum  likelihood 
estimators  (MLEs)  and  descriptions  of  Monte  Carlo  simulations  follow  this. 

By  conparing  performance  results  to  those  of  the  Conplex  Ambiguity  Function 
(CAF),  it  is  shown  that  the  new  MLE(s)  outperform  the  CAF  in  nearly  all  cases.  Also,  by 
conparison  to  the  Cramer-Rao  lower  bound,  the  new  MLE(s)  provide  a  hint  of  the 
ultimate  power  of  cyclostationary  techniques  over  stationary  techniques,  exceeding  this 
stationary  theoretical  lower  bound  on  performance  in  many  cases.  Finally,  the  new 
MLE(s)  perform  comparably  in  the  measurement  of  TDOA  to  those  cyclostationary 
algorithms  previously  derived,  when  adjustments  are  made  for  the  different  signal 
environments  used  in  the  heritage  work. 


VI 


TABLE  OF  CONTENTS 


l.  INTRODUCTION . 1 

A.  MOTIVATIONS  FOR  EMITTER  LOCATION . 1 

1.  Civil  and  Military  Application . 1 

2.  The  Changing  Signal  Environment . 1 

B.  TECHNICAL  PROBLEM  DESCRIPTION . 4 

1.  Problem  of  Interest . 4 

2.  Signal  Models . 5 

C  OBJECTIVES  OF  THE  RESEARCH . 10 

1.  Maximum-Likelihood  Cyclostationary  TDOA  and  FDOA  Algorithms . 10 

2.  Lower  Bound  Selection  and  Computation . 10 

3.  Monte  Carlo  Testing  and  Comparison . 10 

D.  TECHNICAL  APPROACH . 11 

1.  Fraction-of-Tin»  Probability  Approach . 1 1 

2.  Algorithm  Development . 12 

E.  DISSERTATION  ORGANIZATION . 12 

1.  Part  I:  Underlying  Theories . 12 

2.  Part  H:  Development  of  New  Algorithms . 13 

3.  Part  ni:  Monte  Carlo  Testing  and  Conclusions . 13 

n.  TRADITIONAL  TDOA  AND  FDOA  METHODS . 15 

A  TRADITIONAL  GENERALIZED  SIGNAL  MODEL . 15 

1,  Model  B  Revisited . 15 

2.  Comments  on  Simplifying  Assumptions . 16 

B.  GENERALIZED  CROSS-CORRELATION  FUNCTION . 17 

1.  Derivation . 17 

2.  Estimation  Variance . 20 

3.  Binary  Phase  Shift  Keying  Example . . . 21 

C  COMPLEX  OR  CROSS  AMBIGUITY  FUNCTION . 25 

1.  Derivation . 25 

2.  Estimation  Variances . 26 

3.  Binary  Phase  Shift  Keying  Example  Extended . 27 

m.  CYCLOSTATIONARY  THEORY . 35 

A  BACKGROUND . 35 

1.  Basic  Concept . 35 

2.  Philosophies  on  the  Estimation  of  Parameters . 37 

B.  CYCUC  CORRELATION  FUNCTIONS . 42 

1.  Continuous-Time  Cyclic  Correlation  Functions . 42 

2.  Discrete-Time  Cyclic  Correlation  Functions . 45 

C  SPECTRAL  CORRELATION  DENSITY  FUNCTIONS . 46 

1.  Continuous-Time  Spectral  Correlation  Density  Functions . 46 

2.  Discrete-Time  Spectral  Correlation  Functions . 49 

D.  DIGITAL  COMPUTATION  OF  THE  SCD . 50 

1.  Introduction . 50 

vii 


2.  Strip  Spectral  Correlation  Analyzer . 51 

IV.  CYCLOSTATIONARY  TDOA  AND  FDOA  ESTIMATION . 57 

A.  CYCLOSTATIONARY  TIME  DIFFERENCE  OF  ARRIVAL . 57 

1.  Cyclostationary  TDOA  Signal  Model . 57 

2.  Current  Cyclostationary  TDOA  Determination  Algorithms . 59 

V.  NEW  MAXMUM-LIBOELIHOOD  ESTIMATORS .  65 

A.  BACKGROUND . 65 

1.  Signal  Model . 65 

2.  Mathematical  Model . 66 

B.  DERIVATION . *67 

C.  ALGORITHMS . 79 

1.  Approximated  2-Dimensional  Maximum  Likelihood  Algorithm  #1 . 79 

2.  Approximated  2-Dimensional  Maximum  Likelihood  Algorithm  #2 . 81 

VI.  MONTE  CARLO  PERFORMANCE  TESTING . 85 

A.  OBJECTIVES . ”."85 

B.  TEST  ENVIRONMENTS . 86 

1.  Introduction . g6 

2.  Wideband  Interference  Environment . 87 

3.  Narrowband  Interference  Environnaent .  89 

C.  TESTING . 

1.  Stability  of  Estimates . 91 

2.  Products  of  Monte  Carlo  Testing . 93 

D.  RESULTS . 97 

1.  Performance  in  Environment  I  -  Wideband  Interference . 97 

2.  Performance  in  Environment  n  -  Narrowband  Interference . 101 

Vn.  OBSERVATIONS  AND  CONCLUSIONS . 107 

A.  COMMENTS  ON  MONTE  CARLO  RESULTS . 107 

1.  Ideal  vs.  Approximation . 107 

2.  New  Cyclostationary  Estimators  and  the  Cross  Ambiguity  Function . 108 

3.  New  Cyclostationary  Estimators  and  the  Cramer-Rao  Bounds . 109 

4.  Normalization  and  Conparison  to  Heritage  Algorithms . 110 

5.  Another  Look  and  Performance  vs.  Collect  Length . Ill 

6.  Performance  as  a  Function  of  Signal  to  Interference  and  Noise  Ratio . 1 17 

7.  Summary . I33 

B.  ORIGINAL  CONTRIBUTIONS . 135 

C.  AREAS  FOR  FURTHER  INVESTIGATION . 136 

1.  Time-varying  Doppler . 136 

2.  Various  Spectral  Correlation  Function  Substitution  Options . 137 

3.  Amplitude  and  Phase  Mismatch . 137 

4.  Modulation  Types . 138 

5.  Other  Signal  Models . 138 

6.  Lower  Bound  Extensions . ^ . 138 

7.  Other  Optimization  Techniques . 138 

viii 


LIST  OF  REFERENCES . 139 

BIBUOGRAPHY . 141 

APPENDIX . 143 

INITIAL  DISTRIBUTION  LIST . 145 


ix 


X 


ACKNOWLEDGEMENTS 


This  journey  started  for  me  in  1994  when  I  took  my  first  digital  logic  course  and 
was  hooked  on  electrical  engineering.  Now  several  years,  one  dog,  one  son,  two  houses, 
two  jobs  and  a  quarter-of-a-miflion  airline  miles  later,  there  are  a  great  many  people 
without  whom  this  would  not  have  been  possible.  Hopefully  I  have  managed  to  include 
them  all  here. 

First  to  my  committee  I  send  a  heartfelt  thanks.  Having  a  nagging,  email-sending, 
phone  message-leaving  Ph.D  candidate  that  is  3000  miles  away  (makes  strangulation 
difficult)  who  tries  to  demand  very  fast  responses  and  tough  scheduling,  could  not  have 
been  easy.  I  thank  you  all  in  general  for  your  tolerance  of  this  unique  and  certainly  at 
times,  finstrating  endeavor. 

More  specifically,  to  Dr.  Robertson,  thanks  for  your  easy  Texas  manners  and  your 
amazing  ability  to  really  teach  the  complicated  mathematics  that  makes  up  modem  digital 
communications  processes.  I  actually  chose  communications  as  my  major  area  because  of 
what  I  learned  from  you  in  aU  those  classes  and  well-spent  hours  in  the  lab. 

To  Dr.  Fonts-  thanks  for  your  he^)  during  my  Written  Screening  Preparation  as 
well  as  your  calm  and  logical  examination  of  me  during  both  my  Oral  Screening  Exam  and 
Oral  Qualifiers.  I  never  ceased  to  learn  a  new  point  of  view  when  preparing  for  an  exam, 
reviewing  problems  or  answering  your  examination  questions.  In  addition,  thanks  for 
always  working  Space  Systems  Engineering  into  what  we  were  doing  -  the  minor  area 
should  never  be  forgotten  and  you  were  sure  that  it  was  not. 

My  gratitude  extends  to  Dr.  Brown  for  giving  more  of  his  time  and  attention  than 
he  probably  thought  would  be  needed  when  he  agreed  to  do  this.  With  you,  I  was  able  to 
learn  not  to  fear  working  from  first  principles  -  a  lesson  I  will  always  remember.  Without 
you,  never  would  I  have  been  able  to  wade  through  the  reams  of  previously  derived 
cyclostationary  algorithms  that  were  the  basis  for  this  work. 

I  am  indebted  to  Dr.  Lott  for  the  many  hours  of  general  counseling  and  the 
seemingly  endless  answers  he  had  to  the  equally  endless  list  of  questions  I  had  for  him.  I 
found  that  it  is  important  to  have  one  person  on  your  committee  that  you  can  always  call 


xi 


(at  any  hour)  and  ask  any  question  no  matter  how  trivial  or  stupid  it  may  seem  —  you  were 
he.  It  is  rare  indeed  when  a  student  gets  the  opportunity  to  interact  with,  on  a  regular 
basis,  someone  who  has  seemingly  done  everything  there  is  to  do  in  electrical  engineering. 
I  know  of  no  other  person  that  has  accomplished  so  much  at  such  a  young  age  and  still  has 
the  time  to  be  such  a  great  mentor.  I  hope  to  continue  to  learn  from  your  wealth  of 
experience  for  many  years. 

Finally,  I  am  truly  grateful  to  Dr.  Loomis.  For  over  five  years  you  have  expertly 
guided  me  through  a  maze  of  academic  achievements  -  all  of  which  I  never  dreamed  I’d 
see.  Phenomenally,  we  have  managed  to  squeeze  out  three  graduate  degrees  in  that  time  - 
all  because  you  always  knew  exactly  what  to  do,  where  I  should  go  and  who  I  needed  to 
see.  You  always  gave  me  an  incredible  amount  of  latitude  -  certainly  enough  rope  with 
which  to  hang  myself  if  I  so  chose.  However,  you  never  let  me  quite  drop  and  I  learned 
an  enormous  amount  about  how  to  think,  how  to  work  and  how  to  keep  all  the  balls  in  the 
air  at  the  same  time.  Thank  you  sir,  for  all  those  years  and  all  that  time.  I  look  forward  to 
working  with  you  in  the  years  to  come. 

Now  to  the  many  people  that  also  played  in^ortant  roles  I’d  also  like  to  express 
my  gratitude.  Thanks  to  Dr.  Charles  Therrien  for  his  assistance  with  discrete  jointly- 
Gaussian  pdfs  and  for  always  having  time  for  me  to  ask  questions  of  a  true  giant  in  the 
field  of  signal  processing  and  random  processes;  to  Cal  Larvick  for  his  help  on  the 
commercial  uses  of  this  research  area  and  for  reading  the  early,  extremely  rough  draft  of 
this  dissertation  and  providing  very  thorough  comments;  to  my  friends  Val  Johnson  and 
Mitch  Darling  (for  some  major  CMA)  and  the  men  and  women  of  ACG  who  encouraged 
me  throughout;  to  the  GeoLITE  program  office  for  epitomizing  the  notion  that  it’s  what 
you  do  with  it  once  you  have  it  that  counts  in  the  end  and  specifically  to  my  boss,  Randi 
Moore  for  her  constant  support  of  this  work. 

Finally,  there  have  been  several  other  people  without  whom  my  aradftmir  career 
would  have  ended  in  December  1996  with  my  departure  from  NPS.  First,  Linda 
Whitaker,  my  first  post-graduate-school  boss,  has  since  taught  me  everything  I  know 
about  acquisition,  perseverance  through  her  personal  example,  some  humility  along  the 
way  and  the  absolute  importance  of  taking  care  of  your  people.  At  times  she  played 
mother,  boss  and  sister  but  always  friend.  Without  her  unquestioning  support,  constant 

xii 


encouragement  and  the  unprecedented  tune  she  gave  me,  this  work  would  never  have 
been  completed.  Thank  you  Linda  -  we  will  surely  work  together  again  someday. 

To  my  great  friend  Gus  Lott  and  his  wonderful  wife  Mary  -  you  are  two 
tremendous  human  beings  that  I  have  cherished  since  the  day  we  met.  You  have  helped 
me  in  many  ways  well  beyond  careers  and  the  pursuit  of  graduate  degrees.  Gus,  you  got 
me  started,  kept  me  going  and  never  gave  up  on  me.  You  taught  me  and  continue  to 
teach  me  something  new  about  life  every  time  I  talk  to  you.  After  all,  you  are  the  first  guy 
I  found  as  weird  as  I  am  when  it  comes  to  learning  things  simply  from  reading  books! 

Your  influence  extends  far  beyond  circuits  and  antennas  -  you  have  truly  helped  shape  my 
career  and  helped  me  to  realize  what  I  want  to  be  when  I  grow  up  someday.  At  times  I 
still  wonder  if  your  steadfast  support  of  me  in  those  early  years  wasn’t  to  you  own 
detriment.  In  the  end  though,  I  guess  you  could  say  we  are  both  better  off  anyway... 
Mary,  you  certainly  didn’t  need  (and  are  not  nearly  old  enough)  to  have  a  30  year-old  kid 
-  but  you  adopted  one  anyway.  Thanks  for  playing  Mom  and  for  being  the  beautiful, 
caring  fiiend  that  you  are.  Your  zest  for  life  and  your  constant  positive  outlook  on 
everything  is  infectious.  Thank  you  both  for  the  great  years  past  and  the  many  better  ones 
to  come. 

Next,  to  the  best  friend  I’ve  had  for  26  of  my  30  years,  what  can  I  say?  There  isn’t 
anything  we’ve  done  within  our  memories  that  the  other  wasn’t  around  for.  This  time 
around  I  have  to  thank  you  for  reminding  me  that  the  good  life  is  often  more  about  red  felt 
pool  tables  and  muddy  mountain  bikes  than  it  is  about  algorithms  and  derivations.  You 
are  and  will  always  be  my  alter  ego  -  remember,  we  still  share  this  brain,  of  that  I  am 
certain. .  .thanks  Rick. 

Lastly,  I  would  be  truly  lost  without  my  amazing  family.  I  am  indebted  to  my 
brother  Qiristopher  and  my  soon-to-be  sister-in-law,  Tina  for  all  the  he^)  they’ve  given 
me.  From  your  tireless  nephew  entertainment  to  taking  care  of  the  house  and  all  its  animal 
inhabitants  while  we’ve  been  gone,  you’ve  been  a  great  help.  I  hope  to  repay  you 
someday  for  all  the  time  you’ve  given. 

To  my  wonderful  mother-in-law,  Carolyn,  I  cannot  possibly  thank  you  for  all 
you’ve  given  me.  The  personal  example  you  set  in  support  of  education  every  day  and  the 
endless  list  of  things  you’ve  done  for  us,  especially  while  I’ve  been  in  school,  are 


xiii 


staggering.  Thank  you  for  your  unconditional  support,  your  prayers  and  most  especially, 
your  daughter.  I  was  truly  blessed  the  day  I  met  Shannon  because  I  was  able  to  gain  not 
only  a  wife  but  also  another  Mom  in  the  process. 

Mom  and  Dad  -  you  have  always,  since  I  can  remember,  stressed  the  importance  of 
education,  of  applying  oneself  to  the  fullest  and  of  pursuing  dreams.  Dad  -  you  taught  me 
how  to  print  zeros  and  sevens  correctly,  the  futility  of  intelligence  without  skilled  written 
and  oral  ejipression  and  that  learning  is  always  superior  to  mere  grades.  Mom  -  you  never 
let  science  and  engineering  completely  consume  me  with  trips  to  see  Carol  Channing  as 
Dolly  and  James  Earl  Jones  as  Othello,  plane  tickets  to  England  and  trips  to  the  Folger. 
You  both  made  many  sacrifices  when  I  was  younger  to  send  me  to  the  right  schools,  to 
support  me  m  the  pool  and  to  keep  me  out  of  trouble.  You  drove  back  and  forth  to 
Annapolis  more  times  than  you  care  to  he^)  me  survive  at  USNA.  You  have  bought  me 
books  (one  by  Dr.  Seuss  in  particular  made  quite  a  difference),  paid  for  tuition,  dog  and 
cat  watched  and  babysat.  Perhaps  most  nobly,  you  have  suffered  through  my  endless 
descriptions  of  how  the  research  is  going,  what’s  left  and  what  doesn’t  work.  I  cannot 
think  of  two  people  who  could  have  gotten  me  here  better  than  you  did.  I  only  hope  that  I 
can  be  half  the  parent  for  Reilly  that  you  have  been  for  me.  Thank  you  for  everything. 

Finally,  I  have  been  blessed  with  a  remarkable  son,  Reilly  and  an  extraordinary 
wife.  Shannon.  To  my  best  buddy,  you  probably  won’t  remember  much  of  this,  but  let  it 
be  known  that  you  possess,  at  the  r^  old  age  of  three,  a  depth  of  understanding, 
reasoning  and  maturity  that  astounds  me  every  day.  You  have  always  understood  when 
Dad  had  work  to  do  and  always  tried  to  find  any  way  to  help.  Though  I  tried  to  limit  the 
impact  on  your  life,  you  certainly  did  not  get  all  of  Dad,  an  of  the  time,  but  that  never 
seemed  to  bother  you.  Couldn’t  have  done  without  you,  this  is  as  much  yours  as  mine  - 
thanks  for  hewing  out  my  friend. 

My  dearest  -  though  I  certainly  have  given  you  enough  reasons  to  do  so,  you  have 
never  doubted  my  ability  for  a  second.  Even  in  my  darkest  moments  you’ve  always 
managed  to  drag  me  up,  dust  me  off  and  shove  me  back  out  into  the  game.  You  have 
tolerated  my  inexplicable  obsession  with  this  insane  exercise  (and  many  others  along  the 
way,  as  weU  as  the  ones  currently  in  the  works)  with  the  grace  and  poise  of  a  queen.  You 
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are  my  one  true  thing,  with  all  that  I  know  that  means  to  you.  I  love  you  very  much  - 
only  72  or  so  to  go.  TTDT!! 
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INTRODUCTION 


A.  MOTIVATIONS  FOR  EMITTER  LOCATION 


1.  Civil  and  Military  Application 

Emitter  location  enjoys  extensive  civil  and  military  use  as  it  contributes  to  a  broad 
range  of  services.  It  finds  use  in  law  enforcement  during  the  tracking  and  recovery  of 
stolen  vehicles.  It  is  vital  to  search  and  rescue  operations  -  seeking  anything  from  sh4)s 
adrift  to  downed  aircraft  and  lost  hikers.  In  navigation  emitter  location  uses  range  from 
the  multitude  of  Global  Positioning  System  (GPS)  based  applications  to  the  LORAN, 
OMEGA  and  TRANSIT  systems  that  proceeded  it.  Finally,  it  has  found  extensive 
employment  in  the  enforcement  of  communications  regulations  against  unlicensed 
transmissions  and  intentional  or  unintentional  interference  of  legitimate  communications 
resources  recently. 

By  far  the  most  well  known  civil  application  is  GPS.  GPS  receivers  use  the 
difference  in  arrival  times  (or  time  difference  of  arrival  -  TDOA)  between  the  signals, 
ideally  fi:om  four  or  more  GPS  satellites,  to  determine  the  user’s  position  on  the  earth  in 
all  three  dimensions  to  within  100  m  routinely. 

In  military  use,  emitter  location  supports  the  three  prongs  of  Information 
Operations:  protect,  exploit  and  attack.  AU  rely  upon  the  knowledge  of  the  enemy’s 
location.  Protection  of  friendly  communications  requires  the  elimination  of  enemy 
jamming  and  interference  equipment.  Locating  such  equipment  is  the  vital  first  step  in  that 
process.  Exploiting  enemy  systems  too  requires  knowledge  of  their  locations.  Finally,  in 
order  to  attack  enemy  ioformation  systems,  location  data  must  be  available. 

2.  The  Changing  Signal  Environment 

Obviously  the  field  of  emitter  location  has  a  multitude  of  applications  and  draws  its 
motivation  from  both  civil  and  military  uses.  It  has  for  many  years  been  the  subject  of 
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intensive  study,  providing  a  profusion  of  information  on  any  of  the  common  techniques 
and  their  subsets  of  hybrids  and  derivatives.  The  various  algorithms  and  implementations 
have  been  developed,  optimized  and  tested  thoroughly,  producing  volumes  of  signifirant 
work  in  these  traditional  areas.  However,  given  the  increasingly  dense  radio  frequency 
(RF)  environments  throughout  the  world,  interference  tolerant  methods  are  becoming 
more  sought  after.  As  the  RF  environmental  density  increases,  the  methods  must  improve 
to  still  meet  the  accuracies  demanded  of  users. 

Fortunately,  a  recent  development  in  signal  processing  offers  the  possibility  of 
counteracting  some  of  the  deleterious  effects  of  these  increasingly  populated  RF 
environments.  By  modeling  signals  as  cyclostationary  and  taking  advantage  of  the  unique 
features  that  these  models  unmask,  newly  developed  signal  processing  algorithms  are  able 
to  be  highly  signal  selective  even  in  the  presence  of  noise  and  interference  normally 
debilitating  to  traditional  techniques.  Emitter  location  methods  benefit  substantially  from 
the  application  of  cyclostationarity,  particularly  in  the  dense  signal  environments  where  the 
signal  of  interest  (SOI)  shares  its  channel  with  many  spectrally,  spatially  and  temporally 
overlapping  signals  not  of  interest  (SNOI). 

Traditional  methods  of  parameter  estimation  for  emitter  location  perform  well  in 
fairly  benign  signal  environments.  That  is,  with  a  SOI  with  a  significant  positive  signal  to 
noise  ratio  (SNR),  in  the  absence  of  jammers  and  other  SNOI(s)  in  the  frequency  band 
around  the  SOI,  the  traditional  estimators  of  angle-of-arrival  (AOA),  frequency- 
difference-of-arrival  (FDOA)  and  time-difference-of-arrival  (TDOA)  are  superb. 

However,  in  a  peacetime  EM  spectrum  increasingly  packed  with  signals,  and  a  batdespace 
spectrum  potentially  flooded  with  SNOI(s),  hostile  jamming  and  weak  SOI(s),  these 
methods  suffer  from  severely  decreased  signal  to  noise  and  interference  ratio  (SNIR)  as 
well  as  little  signal  selectivity  capability. 

Clearly  for  the  future,  the  move  is  to  ever  higher  portions  of  the  frequency 
spectrum,  heavy  frequency  reuse,  interference  problems  currently  within  the  existing 
spectrum  and  the  periodic  reallocation  of  frequencies  from  one  sector  to  another  (e.g. 
military/  government  to  civil/commercial  and  vise- versa).  All  of  these  frequency 
allocations  are  worked  within  the  bounds  of  the  Federal  Communications  Commission 


2 


(FCC)  for  U.S.  domestic  matters  and  via  the  National  Telecommunications  and 
Information  Administration  to  United  Nations  sanctioned  bodies  for  international  issues, 
the  International  Telecommunications  Union.  Additionally,  the  Worldwide  Administrative 
Radio  Committee  (WARC)  meets  in  Geneva,  Switzerland  every  two  years  to  debate  the 
regional  and  global  application  of  all  portions  of  the  RF  spectrum.  It  is  at  these  sessions 
that  many  of  the  key  decisions  concerning  the  use  of  the  ever-smaller  available  portions  of 
the  RF  spectrum  are  made. 

Additionally,  because  the  demand  has  far  outpaced  the  allocations  in  many  venues 
(U.S.,  Europe  and  Japan),  commercial  wireless  communications  vendors  have  moved  to 
advanced  techniques  for  interference  mitigation.  Taking  advantage  of  both  hardwired 
digital  signal  processing  techniques  and  advanced  software,  they  struggle  to  produce 
handsets,  base  stations  and  systems  capable  of  meeting  demand  and  maintaining  quality  of 
service. 

Further,  by  U.S.  and  international  regulation,  an  increasing  awareness  of  the 
benefits  of  mandatory  emergency-distress-calling  emitter  location  has  placed  new  burdens 
of  mplementation  on  the  service  providers.  In  this  case,  if  a  cellular  telephone  caller  in 
the  United  States  makes  a  9-1-1  call  in  an  emergency,  the  service  provider  within  a  very 
short  time,  is  asked  to  provide  location  information  to  emergency  service  crews  to  find 
this  distressed  caller.  Similar  mandates  are  being  imposed  on  international  cellular  (e.g. 
Global  Service  Mobile  (GSM),  IS95  (CDMA  standard),  etc.)  and  also  on  satellite  mobile 
handsets  (e.g.  Iridium,  Globalstar  etc.)[2]  [3].  The  techniques  inq>lied  in  this  dissertation 
may  have  utility  to  commercial  wireless  providers  in  inqplementing  local,  regional,  or 
global  positioning  schemes  within  their  systems. 

Finally,  other  possible  civil,  commercial  and  military  applications  for  these 
techniques  include:  (1)  Cellular/Personal  Communications  System  (PCS)  interference 
mitigation  rejection  -  caller  and  tower  location  information  can  be  used  to  control  smart 
(dynamic)  antennas  on  both  the  tower  and  the  mobile  subscriber  to  block  out  interfering 
signal  sources,  enabling  greater  geo-ffequency  reuse.  (2)  Transportation  Management 
Applications  -  some  systems  such  as  Omnitracs  by  Qualcomm  have  alr^dy  implemented 
within  which  mobile  units  report  out  their  location  to  a  mother  node  by  query  to  GPS 
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receivers.  If  the  inched  techniques  here  were  implemented,  perhaps  the  GPS  component 
would  not  be  required  and  the  communication  system  itself  could  determine  relative 
location  of  the  mobile  units,  eliminating  the  need  for  anything  but  a  radio  on  the  mobile 
units  and  smart  software  in  the  base  nodes.  (3)  Within  the  U.S.  Department  of  Defisnse, 
Blue  (friendly)  force  location  services  could  be  managed  using  automated  emitter  location 
techniques.  (4)  Finally,  improved  tracking  of  attached,  small  beacons  on  the  animals  under 
intensive  study  could  support  the  tracking  of  migratory  animals  such  as  whales,  birds  and 
polar  bears.  This  is  done  today,  in  a  crude  sense,  using  the  U.S.  and  French  Argos 
satellite  but  perhaps  could  be  more  easily  implemented  using  simple  satellite 
communications  devices  and  service  from  two  or  more  satellites[l]. 

B.  TECHNICAL  PROBLEM  DESCRIPTION 


1.  Problem  of  Interest 

The  siir^lified  technical  problem  addressed  by  this  research  is  that  of  estimating 
TDOA  and  FDOA  for  SOIs  in  a  co-channel  interference  (CCI)  environment  for  the 
purpose  of  emitter  location.  The  relative  difference  between  correlated  signal  conq)onents 
both  in  the  time  domain  as  well  as  the  frequency  domain  at  two  or  more  spatially- 
separated  sensors  makes  possible  the  determination  of  the  direction-of-arrival  (DOA)  of 
the  signal  Those  sensors  may  be  located  on  the  same  platform  or  disparate  platforms. 
Sensors  on  the  same  platform  tend  to  employ  phase-difference  techniques  to  calculate  the 
angle  of  arrival  of  the  SOI  -  a  technique  related  to  the  DOA  solution  from  above. 
Nonetheless,  with  multiple  DOA(s)  from  any  method,  it  is  then  possible  to  determine  the 
location  of  the  emitter. 

The  challenge  in  this  signal  in  interference  environment  and  the  focus  of  this 
research  lies  in  correctly  selecting  the  SOI  among  many  SNOI  and  in  correctly 
determining  the  TDOA  and  FDOA  of  the  SOI,  in  collection  environments  where  negative 
SNIR  is  common.  Traditional  methods  are  susceptible  to  spatial  spectral  and  temporal 
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overlap  rendering  the  task  of  geolocating  a  weak  emitter  among  many  strong  emitters 
difScult  and  at  times  imprecise.  Thus,  a  method  that  is  both  highly  signal-selective  and 
operates  with  high  fidelity  in  poor  SINR  scenarios  is  required. 

The  estimators  derived  herein  gain  their  heritage  from  a  revolutionary  approach 
begun  by  Dr.  William  A.  Gardner  and  his  colleagues  (see  reference  list  and  bibliography 
further  reading).  They  exploit  a  statistical  property  in  the  man-made  signals  known  as 
cyclostationarity.  This  cyclostationarity  is  conspicuous  in  both  the  time  domain  and  the 
frequency  domain  manifesting  itself  as  periodic  components  in  the  mean  and 
autocorrelation  function  (as  well  as  higher  order  moments)  in  the  time  domain,  as  well  as 
power  spectral  density  function  (again  as  well  as  higher  order  moments)  in  the  frequency- 
domain.  Often  the  cyclostationary  statistical  characteristics  of  each  signal  in  a  bandwidth 
are  unique,  providing  an  identifying  signature  that  enables  outstanding  signal-selectivity. 

In  addition,  the  cyclostationary  characteristics  used  for  such  a  high  degree  of  signal 
selectivity  lend  themselves  to  extracting  very  low  SNIR  signals  from  a  highly  corrupt 
environment.  To  that  end,  cyclostationary  signal  processing  techniques  have  been 
developed  for  numerous  purposes  including  weak-signal  detection,  signal  separation  in  co¬ 
channel  interference,  array  processing,  spatial  and  temporal  filtering,  modulation 
classification  and  of  course  TDOA  and  FDOA  estimation. 

2.  Signal  Models 

This  effort  revolves  around  a  series  of  increasingly  realistic  and  complicated  signal 
models  for  the  location  of  an  EM  emitter  through  passive  means  based  upon  measuren^nt 
of  both  the  TDOA  and  FDOA  of  the  SOI.  In  Spooner  and  Gardner  [4],  the  signal  models 
begin  with  a  single  two  receiver,  single  SOI  in  Gaussian  noise  model  This  model 
represents  the  suiq)lified  situation  that  exists  when  the  SOI,  radiating  from  a  source 
remote  to  both  receivers,  impinges  upon  two  spatially  separated  receivers  with  no  relative 
motion  in  the  system.  In  addition,  the  SOI  is  received  in  the  presence  of  Gaussian  noise 
only.  The  observables  for  this  model,  denoted  Model  A,  are 
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x(t)=s(t)+nx(t) 


(1-1) 


y(t)=s(t-D)+ny(t)  (1-2) 

Note  that  s(t)  represents  the  single  SOI  and  »x(t)  3nd  ny(t)  are  zero-mean  Gaussian 
noise  terms.  The  parameter  of  interest  in  this  model  is  D  the  TDOA.  While  this  model 
includes  some  of  the  most  important  aspects  of  the  situation  and  has  been  used  extensively 
to  smq)lify  the  mathematics  of  past  work  in  TDOA,  it  makes  some  critical  simplifying 
assun^tions  that  reduce  its  applicability  to  real-world  situations. 

First,  Model  A  makes  no  allowance  for  anq)litude  and  phase  mismatch  at  the 
receivers.  Complete  coherence  between  receiver  platforms,  a  situation  not  always 
encountered,  is  presupposed  without  regard  to  the  possible  degradation  in  performance 
the  mismatch  can  have  on  TDOA  and  FDOA  estimation.  Second,  this  model  allows  for 
no  relative  motion  in  the  collection  space,  discounting  or  assuming  a  lack  of  consequential 
Doppler  at  each  receiver.  Should  there  be  relative  motion  in  the  collection  space,  the 
frequency-shifted  versions  of  the  SOI  likely  will  not  correlate  with  any  significance  using 
traditional  correlation  based  methods  (depending  primarily  on  the  SOI  characteristics  and 
the  induced  Doppler).  Thus,  while  this  omission  is  perfectly  reasonable  in  known  static 
situations,  it  fails  to  generalize  the  model  for  the  collection  space  to  include  a  broader 
range  of  possible  scenarios.  Thirdly,  the  assumption  of  a  single  SOI  alone  in  Gaussian 
noise  is  only  applicable  in  a  limited  number  of  situations.  Either  indeed  one  signal  is 
detectable  amidst  Gaussian  noise  (highly  unlikely  other  than  in  the  laboratory)  or  the  SOI 
has  such  high  SNR  that  it  can  safely  be  considered  dominant  enough  in  the  environment  to 
validate  the  assumption. 

Finally,  in  most  previous  work  (see  [5],  [6]  and  [7]  for  instance),  the  noise  terms 
are  assumed  independent,  identically  distributed,  Gaussian  distributions  and  uncorrelated 
with  the  SOI  as  welL  While  simplifying  the  mathematics  tremendously,  this  assumption 
can  lead  to  unaccountable  errors  in  the  final  analysis  of  the  performance  of  TDOA  and 
FDOA  estimators.  This  unaccountable  error  is  the  result  of  contributions  to  the  error 
budget  made  by  the  cross-terms  in  the  covariance  matrix  that  are  directly  attributable  to 
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the  correlation  among  the  three  variables.  Indeed,  there  exists  many  situations  where  the 
interdependence  of  the  SOI  and  noise  terms  at  each  receiver  cause  them  to  be  correlated 
but  to  such  a  miniscule  degree  that  the  assumption  remains  essentially  valid.  However, 
depending  on  the  relative  distances  of  the  emitter  and  receivers,  the  atmospheric 
conditions,  frequency  of  transmission,  instrumentation  and  innumerable  other  factors,  the 
SOI  and  associated  noise  terms  at  each  receiver  are  correlated  to  a  signiScant  degree 
(despite  the  normalcy  of  the  noise  probability  density),  thus  requiring  special  consideration 
and  invalidating  this  assumption.  These  four  shortcomings  will  be  addressed  by  four 
additional  signal  models,  each  increasingly  complicated  but  more  realistic  than  the  last. 
They  will  be  designated  Models  B,  C,  D  and  E  and  referred  to  throughout  this  work. 

The  four  additional  signal  models  for  the  TDOA  and  FDOA  emitter  location 
problem  described  above  begin  with  a  model  to  accomodate  the  amplitude  and  phase 
mismatch  found  between  receivers.  Denoted  Model  B,  its  observables 


x:(f)  =  ^(0+Wx(0 

(1-3) 

y(t)  =  ArS(t  -  D) + Tiyii) 

(1-4) 

include  the  conq)lex  constant  Ar  which  represents  the  relative  amplitude  and  phase 
mismatch  between  receivers.  Of  course,  this  mismatch  could  be  represented  by  a  general 
complex  exponential,  Ae^^,  but  for  the  sake  of  simplicity  in  notation,  the  complex 
constant  serves  the  same  purpose  in  a  more  fundamental  fashion. 

Unfortunately,  while  dealing  well  with  the  receiver  mismatch.  Model  B  does  not 
address  the  issues  of  Doppler,  multiple  interfering  signals  or  correlated  SOI  and  noise.  In 
order  to  introduce  Doppler  to  the  collection  space.  Model  C  must  be  used.  While 
maintaining  the  TDOA,  amplitude  and  phase  term  and  Gaussian  noise  terms,  it  includes  a 
Doppler  shift  for  the  second  receiver  relative  to  the  first  (an  FDOA).  The  observables  for 
Model  C  are 
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x(t)=s(t)+nxit) 


(1-5) 


yit)  =  ArS(t-D)e^^^^^  +ny(t)  (1-6) 

where  the  addition  of  fj  represents  the  FDOA  in  the  SOI. 

The  fourth  model  takes  into  account  signals  in  the  collection  space  that  are  not  the 
SOI  but  yet  appear  at  both  receivers.  These  dual-receiver  SNOI  can  severely  interfere 
with  both  TDOA  and  FDOA  estimation  (see  next  Chapter)  and  are  commonly 
encountered  in  the  signal  environments  of  today.  Because  they  appear  at  both  receivers, 
they  contribute  unwanted  correlations  between  the  data  sets  of  the  receivers,  corrupting 
the  measurement  of  TDOA  and  FDOA  for  the  SOI.  The  observables  for  this  model. 

Model  D,  become  very  complicated  and  nearly  intractable  for  theoretical  work. 
Nonetheless,  it  is  inportant  to  realize  them.  They  are 


/ 

=  (1-7) 

i=l 

y(t)  =  ArS(t-D)e^^^^^  -H  +n^(?)  (1-8) 

/=! 

where  the  summation  represents  an  unknown  number  of  SNOI,  7,  each  with  their  own 
TDOA,  D,  and  FDOA,  fd,. 

One  final  signal  model  improves  upon  the  realism  of  Model  D  one  step  further. 

This  model.  Model  E,  accounts  for  SNOI  in  the  collection-space  that  appear  at  one  but 
not  the  other  receiver.  In  fact,  many  collection  space  situations  can  exist  where  the 
contributions  fi-om  these  signals  to  the  analysis  domains  results  in  additional  error  in  the 
estimation  of  the  parameters  of  the  SOI,  further  com()licating  the  analysis.  Model  E  can 
be  represented  as 
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(1-9) 


I  K 

40  =  j(f) ^ 

M  Jfc=l 


y{t)^ArS{t-Dy^^^^ (1-10) 

i=l  l=\ 

where  the  addition  of  an  unknown  number,  K  and  L,  of  single-receiver  SNOI,  40  and 
5i(r),  are  added  to  the  received  signals  x(f)  and  yit)  respectively.  In  addition,  note  that 
while  the  single-receiver  SNOI  have  neither  a  TDOA  nor  an  FDOA. 

In  the  research  presented  here,  all  five  models  will  be  used  either  as  illustrations  in 
the  presentation  of  underlying  theory,  in  the  derivation  of  new  TDOA  and  FDOA 
estimators,  in  the  extension  of  lower  bounds  or  in  the  description  of  the  collection  space 
used  to  test  the  new  algorithms  and  lower  bound  extensions.  While  Models  B  and  C  are 
used  extensively  in  the  theoretical  work  contained  herein.  Models  D  and  E  more 
adequately  represent  the  signal  environment  of  real-world  situations.  Thus,  they  are  used 
as  guidelines  in  the  synthesis  of  test  signals  and  to  represent  real-world  data.  Awareness 
of  the  differences  in  the  signal  models  allows  for  intelligent,  qualitative  comments  to  be 
made  regarding  the  effects  of  using  estimators  derived  for  use  in  a  collection  environment 
characterized  by  Model  C  but  used  in  an  environment  characterized  by  Models  D  or  E  for 
instance.  While  not  quantitative,  the  conclusions  are  valuable  and  may  lead  to  further 
work  that  will  aid  in  optimizing  algorithms  for  use  in  environments  like  those  represented 
by  Model  E. 
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C.  OBJECTIVES  OF  THE  RESEARCH 


1.  Maximum>Likelihood  Cydostationaiy  TDOA  and  FDOA  Algorithms 

The  primary  objective  presented  here  is  the  development  of  two  new 
cyclostationary  TDOA  and  FDOA  estimator  through  the  maximization  of  a  likelihood 
function.  Through  a  new  derivation  with  a  unique  signal  model  and  the  multitude  of 
options  that  are  available  with  the  selection  of  valid  mathematical  substitutions,  one  ideal 
and  two  approximated  maximum-likelihood  estimators  have  emerged. 

2.  Lower  Bound  Selection  and  Computation 

A  secondary  objective  involves  the  survey,  selection  and  possible  modification  of 
an  appropriate  theoretical  lower  bound  (LB)  on  the  estimates  for  TDOA  and  FDOA 
obtained  from  cyclostationary  estimators.  Ideally,  following  a  complete  survey  and 
evaluation  of  related  LBs  from  past  work  in  emitter  location,  the  application  of  an 
appropriate  1^  is  possible.  Such  a  LB  must  be  fully  congruent  with  the  assumptions  and 
basic  theory  of  the  cyclo  stationary  signal  models  and  must  be  easily  calculable. 

Unfortunately,  no  such  LB  exists  currently  despite  of  the  recent  work  by  Schell  on 
the  extension  of  Whittle’s  Theorem  to  the  cyclostationary  case  [8].  Thus,  while  no 
modifications  to  existing  LBs  are  presented  here,  a  comparison  to  the  stationary  Cramer- 
Rao  bound  (CRB)  developed  for  both  TDOA  in  [5]  and  [6]  and  joint  TDOA-FDOA  in 
[7],  serves  to  satisfy  this  secondary  objective. 

3.  Monte  Carlo  Testing  and  Comparison 

The  tertiary  objective  involves  the  thorough  simulation  of  the  estimators’ 
performance  in  the  presence  of  various  challenges  to  traditional  estimators.  Using  a 
Monte  Carlo  method,  designed  both  for  ease  of  comparison  to  previous  algorithm  testing 
and  for  the  assessment  of  specific  signal-selectivity  properties,  the  performance  of  the  new 
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ML  algorithms  is  compared  to  both  the  calculated  theoretical  and  the  performance  of 
previously  derived  estimators.  Performance  conclusions  with  respect  to  several  key  real- 
world  parameters  are  strongly  supported  by  this  testing  and  are  commented  upon  in  the 
final  section. 

D.  TECHNICAL  APPROACH 


1.  Fraction-of-Time  Probability  Approach 

The  theoretical  basis  for  the  models  used  in  the  estimation  of  the  parameters  of 
TDOA  and  FDOA  is  pivotal  to  the  development  of  the  ML  estimators.  Central  to  that 
theory  are  the  two  different  schools  of  thought  related  to  the  estimation  of  parameters 
such  as  TDOA  and  FDOA  fi-om  finite  time-series,  the  traditional  stochastic  approach  and 
the  ffaction-of-time  (FOT)  probability  approach. 

The  traditional  stochastic  process  approach  can  be  traced  through  many  years  of 
use  in  statistical  problems,  particularly  statistical  signal  processing.  This  tradition  makes  it 
a  popular  and  comfortable  choice  for  the  theoretical  underpinnings  of  the  estimation  of 
parameters  process.  However,  it  suffers  m  practical  application  to  actual  finite  time-series 
signal  data  and  fi'om  abstractions  like  the  expectation  operator  that  are  counter-intuitive. 

On  the  other  hand,  the  FOT  probability  approach  defined  by  Gardner  in  [9]  and 
refined  by  Brown  in  [10]  lends  itself  well  to  both  the  analytical  and  empirical  analysis  of 
man-made  signals  represented  by  finite  time-series.  The  FOT  probability  theory 
overcomes  many  of  the  shortcomings  of  the  stochastic  approach  and  avoids  the 
troublesome  abstractions  in  the  process.  For  these  reasons  and  more  (see  Chapter  III),  the 
FOT  theory  serves  as  the  theoretical  basis  for  the  development  of  both  ML  estimators  in 
this  work. 
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2.  Algorithm  Development 


In  previous  work  dating  to  1992,  Gardner  and  his  associates  have  explored  the 
benefits  of  cyclostationary  techniques  in  the  estimation  of  TDOA  and  to  a  lesser  extent, 
FDOA  in  corrupt  signal  environments.  In  [4],  Spooner  and  Gardner  began  the 
development  and  evaluation  of  an  ML  estimator  for  TDOA  based  upon  a  signal  model  that 
included  a  weak  cyclostationary  signal  in  white  Gaussian  noise  with  arbitrary  receiver  gain 
and  phase  mismatch.  In  the  work  presented  here,  an  ML  approach  is  again  pursued  but 
with  a  more  complicated  and  realistic  signal  model,  a  greatly  differing  mathematical  model 
and  several  differing  assun^tions  guiding  the  derivation.  From  that  approach  emerges  an 
ideal  algorithm  optimized  for  both  TDOA  and  FDOA  estimation. 

In  addition,  two  additional  algorithms  are  developed  based  upon  the  investigation 
of  the  individual  terms  of  the  ideal  ML  algorithm.  Making  simple  but  effective 
substitutions  yields  the  first  approximation.  Eliminating  some  weak  terms  to  improve 
confutation  speed  produces  the  second  approximation  and  aids  in  an  initial  evaluation  as 
to  the  feasibility  of  sinf  lifying  the  algorithm  further  to  benefit  implementation  with  little 
effect  on  performance.  Details  of  the  algorithms’  development  appear  in  Chapter  V. 

E.  DISSERTATION  ORGANIZATION 


1.  Part  I:  Underiying  Theories 

In  order  to  appreciate  better,  the  advantages  that  cyclostationary  theory  has  to 
offer  the  emitter  location  problem,  traditional  TDOA  and  FDOA  methods  and  their 
applications  are  briefly  addressed  in  Chapter  n  with  an  enf  basis  on  their  breakdown  in 
highly  corrupt  environments  using  a  Binary  Phase  Shift  Keying  (BPSK)  signal  in  various 
environments  as  a  recurring  example.  Chapter  HI  traces  the  development  of 
cyclostationary  theory  from  the  traditional  auto  and  cross-correlation  fimctions  and  the 
power  spectral  density  function  to  the  cyclic  auto  and  cross-correlation  functions  and  the 
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spectral  correlation  density  function.  This  chapter  serves  as  an  extensive  introduction  to 
the  signal-selective  capabilities  of  cyclostationary  processing  techniques  and  provides  an 
excellent  transition  to  cyclostationary  emitter  location  techniques.  Chapter  IV  introduces 
the  abbreviated  derivation  and  discussion  of  some  current  cyclostationary  TDOA 
techniques. 


2.  Part  n:  Development  of  New  Algorithms 

Chapter  V  offers  the  derivation  of  one  ideal  and  two  approximated  cyclostationary 
joint  TDOA/FDOA  algorithms.  Each  is  derived  in  exhaustive  detail  to  permit  close 
examination  of  the  models  used,  assumptions  made  and  substitutions  used. 

3.  Part  ni:  Monte  Cario  Testing  and  Condusions 

Chapter  VI  tests  the  algorithms’  performance  with  simulated  data  in  a  Monte 
Carlo  fashion,  reporting  the  results  in  terms  of  the  Root  Mean  Square  Error  of  the  TDOA 
and  FDOA  estimates  for  a  variety  of  signal  environments.  In  Chapter  Vn,  by  comparing 
the  results  to  those  obtained  for  the  traditional  Cross  Ambiguity  Function  and  those 
predicted  by  the  stationary  CRB,  the  performance  of  the  new  estimators  relative  to  that  of 
traditional  techniques  is  presented.  Finally,  a  comparison  to  results  from  the  testing  of 
previously  derived  cyclostationary  estimators  precedes  conclusions  and  suggestions  for 
further  research. 
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n.  TRADITIONAL  TDOA  AND  FDOA  METHODS 


A.  TRADITIONAL  GENERALIZED  SIGNAL  MODEL 


1.  Model  B  Revisited 

For  the  development  of  the  two  traditional  noethods  of  TDOA  determination, 
assume  a  stationary  signal  from  a  remote  emitter  impinges  upon  two  spatially  separated 
antenna  elements.  This  situation  may  be  represented  by  signal  Model  B  and  is  repeated 
here  for  consistency  as 

A:(0=5(0+n;c(0  (2-1) 


y(t)=ArS(t-'D)+ny(t)  (2-2) 


Recall  nJif)  and  n/r)  are  assumed  to  contain  only  AWGN  with  no  strong  SNOI(s)  in  the 
frequency  band  of  interest  for  the  SOI,  j(0;  Ar  is  the  relative  magnitude  and  phase 
mismatch  between  the  receivers;  D  is  the  TDOA  of  the  SOI. 

It  follows  that  the  auto-correlation  and  cross-correlation  functions  are  given  by 


(2-3) 

(2-4) 

(2-5) 

and  the  respective  spectral  density  functions  are 
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5,(/)=W)+5,,(/) 


(2-6) 


s,(/)=|^fs,(/)+s^(/) 

(2-7) 

Sy,(f)=ArSs(J)-e-i^’^+Sn^(f) 

(2-8) 

These  relationships  contain  the  parameter  of  interest,  specifically,  D,  the  TDOA. 
The  next  two  sections  illustrate  two  traditional  methods  used  to  extract  D  from  these 
equations.  The  first,  the  Generalized  Cross  Correlation  (GCC)  function,  can  be  used  in 
the  absence  of  relative  motion  between  the  two  receivers  and  the  emitter.  The  second,  the 
Con^)lex  Ambiguity  Function  (CAF),  can  be  used  to  estimate  FDOA  and  TDOA 
simultaneously  as  is  necessary  when  relative  motion  between  the  receivers  and  emitter 
exists  causing  the  use  of  Model  C  to  accommodate  the  Doppler. 

2.  Comments  on  SimpU^ng  Assumptions 

As  alluded  to  in  Chapter  I,  of  key  importance  to  the  development  illustrated  here  is 
the  statistical  independence  of  s(f),  nx(t)  and  ny(t)  and  the  absence  of  in-band  interference. 
While  it  may  be  argued  that  the  signals  as  well  as  noise  measured  at  the  two  receivers  can 
be  correlated  to  some  extent,  that  scenario  greatly  complicates  this  case  and  is  beyond  the 
scope  of  this  discussion.  The  derivations  that  follow  are  drawn  from  [5],  [6],  which  treat 
the  GCC  more  specifically  and  [7]  and  [1 1]  which  address  the  CAF  more  so  than  the 
GCC.  All  assume  the  three  components  to  be  uncorrelated.  For  an  initial  treatment  of 
correlated  noise  samples  at  the  two  receivers  which  produce  non-zero  off-diagonal  terms 
in  the  covariance  matrix  of  the  TDOA  estimate,  see  [12]. 
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B.  GENERALIZED  CROSS-CORRELATION  FUNCTION 


1.  Derivation 

Applying  (2-5)  in  the  absence  of  SNOI  in  the  SOI  band  of  interest  and  assuming 
statistically  independent  s(t),  n/t)  and  ny(t)  yields 


Ry^(T)  =  Ar‘RsiT-D)  +  Rn^ny(.'^)  (2-9) 

This  function  will  peak  at  t=  D,  the  TDOA  between  receivers  1  and  2.  Because  n;^t)  and 
ny(t)  contain  merely  AWGN,  their  cross-correlation  term  in  (2-9)  above  reduces  the  SNR 
of  the  measurement  but  does  not  add  interference  in  the  form  of  SNOI(s).  Further 


analysis  using  (2-6)  -  (2-8)  in  a  similar  manner  reveals 

Sx(f)  =  Ss(f)+Sn^(f)  (2-10) 

Sv(/)  =  KfSs(/)+S,^(/)  (2-11) 

Syx(J)  =  Ar-Sx{f)-e-J'^+S„^„^(f)  (2-12) 

To  extract  D  from  (2-12),  first  note 

W)=-^°  /o-f^/^/o+f  (2-13) 

=  0  otherwise 


which  assumes  the  SOI  effectively  occupies  a  finite  bandwidth  B  around  the  carrier  or 
intermediate  frequency /o.  If  sufficiently  high  SNR  is  assumed,  the  noise  terms  in  (2-10)  - 
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(2-12)  may  be  disregarded.  Thus,  the  ratio  of  the  spectra  using  (2-10)  with  (2-13) 
becomes, 


i  -jlnfD  r  f  f 

SAf)  '  2  ^  2 


(2-14) 


and  taking  the  inverse  Fourier  transform  of  this  ratio  gives 


K(.t)=  J 


/•-y2 


SAf) 


(2-15) 


which  peaks  at  t=  D.  The  noise  terms  in  (2-14)  are  dropped  in  (2-15)  because  they  are 
identically  zero  due  to  the  statistical  independence  of  the  signals  and  the  noise  assumed  in 
this  modeL  Equation  (2-15)  in  turn  may  be  rewritten  as 

fo*% 

KW=  \Wif)Sy^<J)e*i^’r'df  (2-16) 


The  weighting  function,  W(f),  is  defined  in  this  instance  as  USjlfi  in  (2-16)  which 
is  the  best  choice  given  no  a  priori  information  about  the  SOI  [5].  In  addition,  it 
distinguishes  this  case  as  the  generalized  cross  correlation  method.  Assigning  a  weighting 
function,  W(f)  =  1  reduces  the  TDOA  determination  to  a  simple  cross-correlation  as  in  (2- 
5).  Given  prior  knowtedge  of  the  noise  and  interference  characteristics  it  is  possible  to 
make  other  choices  for  W(f)  including  the  Roth  impulse,  SCOT  and  PHAT  among  many 
others  which  are  designed  to  reduce  specific  noise  and  interference  problems  [5].  A  more 
comprehensive  list  of  weighting  functions  and  associated  uses  appears  in  Table  2-1. 
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It  is  clear  that,  in  the  absence  of  significant  noise,  co-channel  interference  and 
relative  motion  between  receivers  and  emitter,  finding  the  peak  in  the  GCC  produces  the 
desired  estimate  of  the  TDOA  of  a  signal  from  the  ratio  of  the  estimates  of  the  power 
spectral  density  function  of  the  signal  at  one  receiver  and  the  cross  power  spectral  density 
function  of  the  two  received  signals.  Figure  2-1  below  illustrates  the  generalized  cross 
correlation. 


Processor  Name  Weighting  Design 

Functions  Goal 


Cross  Correlation 

1 

None 

Roth  Impulse  Response 

A(/) 

Stq>press  noise  in  receiver  #1  for 

the  frequency  region 

surrounding  tbe  SOI 

SCOT 

/j 

Suppress  noise  in  both  recdvers 

for  the  frequency  region 

surrounding  the  SOI 

^xx(f)^yy(f) 

PHAT 

Developed  purely  as  an  ad  hoc 

technique 

Eckart 

Maximizes  the  difference  in  the 

mean  correlator  ouq>ut  for  a 

signal  present  versus  just  noise 

present  (deflection) 

Table  2-1  -  Candidate  weighting  functions,  from  [5] 


I 


19 


Figure  2-1  -  GCC  block  diagram 


2.  Estimation  Variance 

Beginning  in  [5],  Knapp  and  Carter  developed  an  expression  for  the  variance  of 
the  TDOA  estimate  resulting  from  the  GCC.  In  [6],  Azaria  and  Hertz  refined  the 
expression  relating  the  variance  of  the  estimate  to  the  integration  time,  signal  bandwidth 
and  equivalent  SNR  at  each  receiver.  Given  those  inputs, 


,2 

£:(d-d) 


3  111 


(2-17) 


A 

where  D  is  the  TDOA  estimate;  D  is  the  actual  TDOA;  T  represents  the  integration  time; 
Bs  the  signal  and  noise  bandwidths;  —  is  defined  as  the  equivalent  SNR  as  follows 

r 
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(2-18) 


1_1 

1 

y  1 

y  2 

Ti  ri 

nr2 

3.  Binary  Phase  Shift  Keying  Example 

To  best  illustrate  the  performance  of  the  GCC  in  both  benign  and  corrupt 
environments,  assume  a  single  SOI  in  AWGN  as  represented  by  signal  Model  B.  A  BPSK 
signal  with  the  characteristics  listed  in  the  Appendix  serves  as  the  SOL  The  TDOA  of  this 
signal  is  clearly  50  samples  as  the  plot  of  the  output  of  the  GCC  function  for  this  case 
(GCC- 1)  illustrates  in  Figure  2-2. 
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Figure  2-2  -  Plot  of  GCC  output,  high  SNR  (GCC-1) 


For  comparison,  change  the  SNR  of  the  BPSK  signal  to  -10  dB.  The  output  of 
the  GCC  for  this  case  (GCC-2)  appears  in  Figure  2-3.  Note  that  in  the  presence  of  severe 
AWGN,  the  GCC  begins  to  produce  errors  in  the  estimate  of  the  TDOA. 
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Thirdly,  the  BPSK  signal  is  restored  to  its  original  SNR  but  now  has  an  FDOA  of 
20  kHz  introduced  to  the  SOI,  measured  at  the  two  receivers  as  in  Model  C.  Figure  2-4 
illustrates  the  results  produced  by  the  GCC  in  this  situation  (GCC-3).  The  shifting  of  the 
frequency  con^onents  of  the  SOI  was  enough  to  produce  errors  in  the  estimation. 


so  100 

Time  Delay  (samples) 


Figure  2-3  -  Plot  of  GCC  output,  SNR  =  -10  (GCC-2) 

As  a  fourth  illustration,  the  environment  of  GCC-2,  low  SNR  and  GCC-3, 
Doppler,  are  combined  and  the  results  displayed  in  Figure  2-5  (GCC-4).  Note  the 
combination  of  poor  SNR  and  Doppler  results  in  extremely  poor  performance  for  the 
GCC.  In  the  next  section,  the  effects  of  Doppler  are  successfully  accounted  for  in  the 
second  of  the  two  illustrative  examples  of  traditional  stationary  methods  of  TDOA  and 
FDOA  determination,  given  sufficient  SNR. 

Finally,  to  further  illustrate  the  breakdown  of  the  GCC  in  modem  signal 
environments,  note  the  effect  of  introducing  a  single  interfering  SNOI  with  significant 
spatial,  temporal  and  frequency  overlap.  The  single  SNOI  is  another  BPSK  signal  with 
characteristics  shown  in  the  Appendix.  The  effects  of  this  corraption  on  the  estimate  of 
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the  TDOA  of  the  SOI  can  be  seen  in  Figure  2-6  (GCC-5).  Note  the  SOI  and  SNOI  have 
different  center  frequencies,  bandwidths,  keying  rates  and  TDOA(s)  which  might 
distinguish  them  from  one  another  quite  well  under  other  circumstances.  However,  in 
correlating  the  outputs  of  the  two  receivers,  the  combination  of  the  partial  overlap 
spatially,  temporally,  and  spectrally  produces  a  TDOA  estimate  that  represents  the  true 
TDOA  of  neither  emitter  at  109  samples.  In  addition,  though  it  is  obvious  that  two  signals 
are  present,  by  the  two  peaks  in  the  plot,  even  assuming  the  second  peak  at  29  samples  is 
correct  would  lead  to  an  incorrect  estimate  for  the  TDOA  of  the  SOI.  As  a  logical 
extension,  it  is  obvious  that  the  addition  of  poor  SNR  (and  consequently  even  worse 
SNIR)  as  well  as  Doppler  for  the  SOI  to  this  last  situation  would  render  the  GCC 
completely  ineffective. 

While  the  CAF  discussed  in  the  next  section  accounts  for  Doppler,  it  too  suffers 
from  poor  SNR  environments.  In  addition,  it  will  become  evident  that  SNOI  make  TDOA 
(and  FDOA)  estimation  difficult  for  both  the  GCC  and  CAF  and  can  be  handled  only 
through  different  techniques  like  the  exploitation  of  cyclostationarity. 


Figure  2-4  -  Plot  of  GCC  output,  high  SNR  and  Doppler  (GCC-3) 
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Figure  2-5  -  Plot  of  GCC  output,  low  SNR  and  Doppler  (GCC-4) 


C.  COMPLEX  OR  CROSS  AMBIGUITY  FUNCTION 


1.  Derivation 

The  coirq)lex  ambiguity  function  (CAP),  also  known  as  the  cross  ambiguity 
function  or  the  radar  function  can  be  interpreted  as  an  extension  of  the  GCC  for  moving 
transmitters  and/or  receivers.  Since  the  GCC  function  relies  upon  the  correlation  of  the 
signals’  characteristics  in  either  the  time  or  frequency  domain.  Naturally,  if  components  of 
one  of  the  received  signals  were  shifted  in  frequency  by  some  amount  relative  to  the  other 
received  signal  (Doppler),  the  correlation  would  surely  be  corrupted  and  may  even  be 
rendered  useless.  This  characteristic  of  the  GCC  was  illustrated  in  Figures  2-6  and  2-7 
above. 

In  fact,  Stein  in  [7],  showed  that  in  order  to  properly  determine  the  TDOA 
between  two  receivers  in  the  presence  of  a  Doppler  difference  at  each  receiver,  the 
spectrum  of  one  of  the  signals  first  must  be  translated  in  frequency  by  an  amount  fa  equal 
to  the  Doppler  difference  (FDOA)  measured  between  the  observers.  In  order  to  show 
this,  a  Doppler  shift  is  introduced  into  Model  B  to  obtain  Model  C.  Equations  (2-1)  and 
(2-2)  can  now  be  written  as 


(2-19) 


y(t)  =  ArS{t-Ty)e^'^'^^^  +ny{t)  (2-20) 


where  the  addition  of  fd  is  the  FDOA  of  the  SOI  as  measured  between  receivers  1  and  2. 

In  place  of  calculating  TDOA  with  the  two  dimensional  GCC  (time  and  magnitude 
of  the  correlation)  from  (2-16),  which  in  the  presence  of  significant  Doppler,  could  peak  at 
a  value  that  does  not  truly  correspond  to  the  TDOA,  it  is  necessary  to  calculate  the  three 
dimensional  (time,  frequency  and  magnitude  of  the  correlation)  CAF  given  by 
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(2-21) 


0 

A  block  diagram  of  the  CAF  operation  in  the  time  domain  appears  in  Figure  2-7. 


Figure  2-7  -  CAF  block  diagram 


The  simultaneous  determination  of  the  TDOA,  D  and  the  Doppler  shift,  U,  causes 
I  A(D,fj)|  to  peak.  At  =  0,  the  CAF  reduces  to  a  GCC  problem  as  outlined  above.  For 

0,  the  CAF  may  be  thought  of  as  a  GCC  performed  after  frequency  shifting  the 
spectrum  of  y(t)  up  or  down  as  necessary  by  an  amount  equal  to  fd 

2.  Estimation  Variances 

The  three-dimensional  width  of  the  correlation  lobe  in  (2-20)  is  directly 
proportional  to  the  accuracy  of  the  estimates  of  TDOA  and  FDOA.  Stein  further  points 
out  in  [7]  that  the  variance  of  the  estimate  for  each  parameter  can  be  related  to  the  noise 
bandwidth,  the  signal  bandwidth,  the  integration  time  and  the  effective  input  signal  to 
noise  ratio  as  follows 
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_  1  1 

<^TOOA 


(2-22) 


3.  Binaiy  Phase  Shift  Keying  Example  Extended 

For  ease  of  comparison,  the  following  five  examples  apply  the  CAF  to  the 
precisely  the  same  situations  represented  by  the  five  GCC  plots  from  the  previous  section. 
(Reference  the  Appendix  for  the  characteristics  of  the  SOI  and  the  SNOI.) 

In  the  first  two  sets  of  CAF  plots.  Figures  2-8  (CAF-1)  and  2-9  (CAF-2)  show 
three  plots  each.  The  first  plot  is  the  three  dimensional  CAF  plot  while  the  remaining  two 
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are  plots  of  the  maximum  values  along  lines  of  time  delay  and  frequency  delay 
respectfully.  Figure  2-8  represents  performance  in  a  high  SNR  environment  with  no 
Doppler  and  Figure  2-9  illustrates  the  degradation  in  a  low  SNR  environment.  Note  that 
because  the  CAF  reduces  to  the  GCC  in  this  case  only  the  three-dimensional  rendering  and 
slight  differences  in  the  plotting  of  the  time  delay  distinguish  these  plots  from  Figures  2-2 
and  2-3. 

Furthering  the  CAF  example,  introduce  the  same  FDOA  of  20  KHz  used  in  the 
GCC-3  and  GCC-4  cases  above.  Here,  in  the  high  SNR  environment,  the  CAF 
successfully  estimates  both  the  TDOA  and  FDOA  in  Figure  2-10  (CAF-3).  However,  the 
low  SNR  environment  deteriorates  the  accuracy  of  both  estimates  which  may  be  seen  in 
Figure  2-11  (CAF-4).  No  single  feature  dominates  the  estimation  of  either  TDOA  or 
FDOA  and  in  fact,  the  maximum  values  atong  each  axis  are  incorrect  estimates  for  the 
actual  TDOA  and  FDOA  of  the  SOI. 

Finally,  the  introduction  of  the  BPSK  SNOI  in  Figure  2-12  has  a  similar  effect  on 
the  CAF  as  it  has  on  the  GCC.  Note  that  the  SNOI  signal  has  been  modified  to  contain  an 
FDOA  of  -1 1.7  kHz  to  he^  in  the  illustration  of  the  effect  of  the  SNOI  on  the  FDOA 
estimate  as  welL  As  was  true  with  case  GCC-5  above,  the  addition  of  a  strong  SNOI 
spectrally,  temporally  and  spatially  close  to  the  SOI  masks  the  peak  of  the  SOI  in  the 
CAF.  The  dominant  feature  in  Figure  2-12  is  actually  attributable  to  the  SNOI  energy, 
making  the  accurate  determination  of  the  TDOA  and  FDOA  of  the  SOI  impossible 
without  some  pre-filtering  like  that  performed  by  the  weighting  functions  of  the  GCC. 

Low  SNR  and  overlapping  SNOI  make  both  the  GCC  and  CAF  much  less 
effective  in  the  successful  estimation  of  the  parameters  of  TDOA  and  FDOA  for  a  single 
SOI  in  the  collection  bandwidth.  The  GCC’s  inability  to  account  for  Doppler  shifts  and 
both  the  methods’  inability  to  distinguish  the  parameters  of  the  SOI  from  overlapping 
SNOI  motivate  the  pursuit  of  methods  that  support  the  accurate  estimation  of  signal 
parameters,  such  as  TDOA  and  FDOA,  in  the  presence  of  poor  SNR  environments  riddled 
with  SNOI. 
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Figure  2-8  -  Plots  of  CAF  output,  SNR  =  3  dB  and  no  Doppler  (CAF-1) 
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Figure  2-1 1  -  Plots  of  CAF  output,  SNR  =  -10  dB  and  Doppler  =  20  kHz  (CAF-4) 
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Figure  2-12  -  Plots  of  CAF  output,  SINR  =  -13  dB,  Doppler  =  20  kHz,  SNOI  (CAF-5) 


The  exploitation  of  cyclostationarity  in  these  signals  makes  possible  the  accurate 
estimation  of  signal  parameters  both  in  poor  SNR  situations  as  well  as  in  the  presence  of 
many  spectrally,  temporally  and  spatially  overlapping  SNOI.  In  the  next  two  chapters  first 
the  general  theory  of  cyclostationarity  and  then  past  work  on  its  specific  application  to 
TDOA  and  FDOA  estimation  is  documented.  In  this,  the  noise  tolerant,  signal  selective 
characteristics  of  these  techniques  will  be  evident. 


m.  CYCLOSTATIONARY  THEORY 


A.  BACKGROUND 


1.  Basic  Concept 

The  concepts  of  cyclostationarity  have  been  examined  in  theory  for  well  over  two 
decades.  Beginning  in  the  late  1960’s,  Dr.  Lewis  Franks  of  the  University  of 
Massachusetts  at  Amhurst  and  Dr.  William  Gardner  of  the  University  of  California  at 
Davis  (UCD)  began  extensive  research  in  the  area  of  cyclostationary  signal  processing. 

Dr.  Gardner  has  since  produced  a  plethora  of  publications  in  all  aspects  of  the  field.  His 
paper  on  the  general  theory  of  cyclostationary  signal  processing,  published  in  the  April 
1991  edition  of  the  IEEE  Signal  Processing  Magazine  [13],  stands  as  the  key  reference  for 
most  aspects  of  the  basic  theory.  As  is  always  the  case  when  one  individual  has 
contributed  so  much  to  an  important  engineering  development,  much  of  the  theory  and 
examples  presented  here  follow  Dr.  Gardner’s  work  closely. 

As  previously  noted,  most  modem  signal  processing  techniques  associated  with 
communications  applications  treat  SOI(s)  as  stationary  random  processes.  However, 
because  most  manmade  signals  are  generated  through  some  repetitive,  periodic  process 
such  as  the  amplitude,  frequency  or  phase  modulation  of  a  sinusoidal  carrier,  the  encoding 
of  data  or  the  encryption  of  a  message,  their  statistics  inevitably  vary  periodically  with 
time.  While  in  many  instances,  receivers  may  successfully  ignore  the  underlying 
periodicity  of  a  manmade  signal,  often  the  detection  of  the  signal  and  the  estimation  of  its 
parameters  is  more  successfully  accomplished  by  modeling  the  signal  as  cyclostationary 
instead  of  stationary. 

As  stated  before,  a  process  with  statistics  that  vary  periodically  with  time  is  termed 
cyclostationary.  A  stationary  random  message  such  as  digital  data  or  analog  voice  is 
modulated,  clocked  or  framed  by  a  periodic  communications  process  through  some 
nonlinear  system  and  a  cyclostationary  signal  results. 
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From  a  strictly  mathematical  point  of  view,  a  cyclostationary  signal  of  order  n  is 
one  that  will  have  additive  sine-wave  components  that  result  in  spectral  lines  in  some  n®* 
order  nonlinear  transformation  of  the  signal.  In  the  case  of  n  =  2,  a  signal  is  said  to  be 
second  order  cyclostationary  if  a  quadratic  transformation  produces  additive  sine-wave 
components  that  generate  spectral  lines.  This  characteristic  may  be  thought  of  as  akin  to  a 
process  bemg  considered  wide  sense  stationary  or  stationary  through  order  2. 

To  continue  the  case  of  n  =  2  more  specifically,  a  signal  x(f)  is  cyclostationary 
with  cycle  frequency  a  if  and  only  if  some  of  its  delay  products  y(t)  =  x(t)x(t  - 1)  produce 
a  spectral  line  at  fi'equency  a.  If  not  all  cycle  frequencies  a  for  which  jc(f)  exhibits 
cyclostationarity  are  harmonics  of  a  single  fundamental  frequency,  thenx:(0  is 
polycyclostationary  [14].  Polycyclostationarity  implies  the  existence  of  more  than  one 
periodicity  in  the  statistics  of  a  signal  This  in  turn  implies  more  than  one  source  of 
periodicity  such  as  the  case  when  clocked-digital-data  phase  modulates  a  sinusoidal  carrier 
in  an  M-ary  phase  shift  keying  scheme  [14]. 

To  further  illustrate  the  concept,  consider  a  signal  jc(f)  that  contains  an  additive 
sine- wave  component  at  frequency  a  and  is  of  the  form 

x(t)=acos(2oa+0)  with  a  ^0  (3-1) 

The  Fourier  coefficient  defined  as 


M^  =  (x(t)e  (3_2) 

where  (•)  denotes  time-averaging,  will  be  non-zero.  Note  that  this  is  of  the  form  of  the 
common  representation  of  the  Power  Spectral  Density  (PSD)  of  x(r)  with  a  spectral  line  at 
-Kx  and  its  mirror  at  -a.  In  simple  terms,  xit)  contains  first-order  periodicity  at  frequency 
a  given  by 


M?f[^(f-CC)+Sif+a)]  (3-3) 
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Now,  considering  contributions  to  x(t)  from  other  sources,  the  total  signal  may  be 
represented  as 


x{t) = a  cosilTtat +0)+ nit)  (3-4) 

where  nit)  can  be  thought  of  as  the  random  energy  outside  that  of  the  SOI.  If  nit)  is 
strong  in  comparison  to  the  SOI  such  that  it  masks  the  sine- wave  components  in  xit)  from 
detection  during  casual  inspection  of  the  waveform,  then  Ac(f)  can  be  thought  of  as 
possessing  hidden  periodicity.  This  hidden  periodicity  can  still  be  exploited  by  using  signal 
processing  techniques  such  as  the  PSD  function  as  noted  above.  However  more  powerful 
cyclostationary  signal  processing  techniques  are  more  sophisticated  and  can  unmask 
periodicities  in  signals  that  are  hidden  even  from  common  traditional  techniques  like  the 
PSD  function. 


2.  niilosophies  on  the  Estimation  of  Parameters 

Central  to  the  results  presented  as  part  of  this  work  is  the  theoretical  basis  for  the 
models  used  to  estimate  the  parameters  of  TDOA  and  FDOA.  At  the  center  of  that  theory 
are  the  two  different  schools  of  thought  related  to  the  estimation  of  those  parameters  from 
time-series.  At  first  glance,  the  approaches  may  appear  to  be  philosophical  in  nature  only. 
However,  the  two  approaches  diverge  considerably  in  their  application  to  practical 
problems.  In  fact,  one  approach  turns  out  to  be  more  intuitive  than  the  other  and  thus 
more  easily  applied  to  processing  actual  signal  data. 

Traditionally,  the  probabilistic  theory  of  statistical  inference  treats  a  single  time- 
series  as  merely  one  of  an  infinite  number  of  sample  paths  in  the  ensemble  of  the  assumed 
cyclostationary  stochastic-process  model  The  ensemble  is  infinite  in  both  numbers  of 
paths  (down  the  ensemble)  and  in  time  (across  the  ensemble).  Both  these  facts,  while 
strictly  mathematical  creations  are  necessary  in  the  stochastic  process  riiodel  and  as 
outlined  in  more  detail  below,  cause  difficulties  in  their  application  to  practical  problem 
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solving.  Reliability  in  parameter  estimations  made  from  just  a  single  time-series  in  this 
ensemble  is  hinged  upon  obtaining  the  same  values  for  estimates  regardless  of  the  sample 
path  chosen.  This  reliance  upon  reliability  down  the  ensemble  in  the  presence  of  merely  a 
single  tum-series  of  actual  data  again  causes  some  practical  problems  in  obtaining 
empirical  evidence  of  such  reliability. 

This  stochastic-process  model  approach  is  empowered  by  great  tradition  and 
familiarity  and  enjoys  the  vast  historic  underpinnings  of  the  extensive  study  of  random 
phenomena.  The  principles  and  underlying  elementary  mathematics  are  well  known  and 
provide  a  strong  basis  for  the  stochastic-process  model  approach.  In  addition,  the 
stochastic-process  model  fits  well  with  many  engineering  and  science  problems  where 
infinite  ensembles  are  perfectly  acceptable  constructs. 

On  the  other  hand,  it  does  not  lend  itself  intuitively,  as  readily,  to  the  development 
of  techniques  for  the  analysis  of  much  of  the  actual  signal  data  processed  regularly.  Its 
abstractions,  like  the  expectation  operator,  and  reliance  upon  infinite  (in  both  dimensions) 
ensembles  tends  to  make  the  method  detached  from  the  practical  study  of  common  signal 
processing  problems  found  in  everyday  work.  For  instance,  many  communications  signals 
are  inherently  single-sample-path  phenomena  (even  in  instances  where  time,  frequency  and 
spatial  diversity  deliberately  appear.  Thus,  these  types  of  problems  derive  no  practical 
value  from  the  stochastic  process  approach  that  requires  assutq)tions  of  infinite  ensembles 
when  but  one  sample  path  is  possible  m  any  one  experiment  without  some  very  high  level 
generalizations  that  in  the  end,  would  remove  most  of  the  fidelity.  In  addition,  in  the  real¬ 
time  or  near  real-time  signal  processing  problems  of  these  single-sample-path  time-series, 
there  lies  a  danger  in  incorrectly  assuming  the  underlying  model  in  the  stochastic 
approach.  Often  these  mistakes  prove  fatal  to  the  solution  in  practice.  Nonetheless,  the 
stochastic-process  model  for  cyclostationarity  offers  some  advantage  in  that  for  theoretical 
work,  the  probabilistic  method  is  often  easier  to  manipulate.  Thus  the  development  of  this 
approach,  while  ultimately  not  used  here,  is  not  without  merit. 

A  random  process  is  defined  as  first  order  cyclostationary  with  period  T  if  its 
probability  density  ftmction  can  be  characterized  by 
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fx{x{t)}=fx[x{t+nT)] 


(3-5) 


where  t  is  real  and  n  is  an  integer.  First  order  periodicity,  of  course,  implies  that  the  mean 
of  the  function  is  periodic  with  period  T 

E\^x{t)  =  E[x{t  -I-  nr)}J  (3-6) 

Given  a  periodic  mean,  the  expected  value  can  be  represented  in  terms  of  a  Fourier  series 
with  the  form 

£{j:(f)}=  (3-7) 

where  the  Fourier  coefficient  mi  can  be  written  as 


Equation  3-8  represents  the  cyclic  mean  and  is  the  manifestation  of  first-order  periodicity 
in  x(t).  It  can  also  be  considered  the  probabilistic  cyclic  mean. 

A  random  process  is  said  to  be  second  order  cyclostationary  if 

£{:c(fi)x(f2)} = fx{x(.h  +nT),x{t2  +nT)]  (3-9) 

where  ti  and  tz  are  real,  n  is  an  integer  and  T  is  the  period.  From  (3-9)  it  can  be  seen  that 
the  autocorrelation  function  of  a  second-order  cyclostationary  process  contains  is  periodic 
with  period  T  and  can  be  written 
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E{x{t)x{t + t)}  =  Rxx{t+ nT,  r) 


(3-10) 


and  is  considered  cyclostationary  in  the  wide-sense.  In  this  instance,  the  Fourier 
coefficients  derived  from  this  periodicity  are 


(3-11) 


and  are  referred  to  collectively  as  the  probabilistic  cyclic  autocorrelation  function. 

In  addition,  the  probabilistic  spectral  correlation  density  (SCD)  function  is  defined 
simpfy  as  the  Fourier  transform  of  (3-11) 

(3-12) 

Note  that  a  more  rigorous  approach  can  be  found  in  Gardner,  [14].  The  simple  derivation 
above  provides  a  basis  for  comparison  with  the  non-probabilistic  approach  preferred  for 
the  work  contained  herein  and  for  which  the  derivation  appears  in  much  greater  detail 
below. 

Given  the  shortcomings  of  the  stochastic  approach,  an  alternative  that  better 
addresses  the  needs  of  single-sample-path  signal  processing  problems  is  needed.  For  this, 
the  time-series  approach  evolved.  Treating  the  time-series  as  a  fraction  of  an  infinitely 
long  single  sanple  path  avoids  many  of  the  abstractions  associated  with  forcing  the  idea  of 
an  ensemble  from  a  stochastic  process  model  onto  measured  signal  samples.  Reliability  of 
estimates  in  this  approach  relies  only  on  the  ability  to  obtain  approximately  the  same 
values  independently  of  the  fraction  of  time  of  the  infinite  time-series  used  to  obtain  the 
estimate.  Inherently,  this  process  model  applies  well  to  the  types  of  problems  described  as 
troublesome  for  the  stochastic-process  approach,  provided  that  the  time-series  approach 
can  account  for  the  necessary  statistical  inference  and  decision  theory  needed  to  develop 
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and  evaluate  signal  processing  algorithms  within  this  framework.  Fortunately  this  has 
proven  true. 

Essentially,  rather  than  assuming  an  infinite  replication  of  the  experiment  to  form 
an  ensemble,  the  time-series  approach  presumes  a  time-invariance  that  allows  the  signal 
data  to  be  extended  in  both  directions  without  bound.  Then,  the  process  model  is  directly 
related  to  the  data  in  the  limit  as  the  data  approaches  infinity.  This  naturally  leads  to  the 
theory  of  fraction-of-time  (FOT)  probability  where  the  statistics  are  related  to  the  fraction 
of  time  the  data  can  be  expected  to  behave  in  a  measured  manner. 

Gardner  in  [15]  proposes  that  the  theory  of  statistical  inference  and  decision  in  the 
FOT  model  is  isomorphic  to  that  of  the  stochastic-process  model.  To  that  end,  the  goal  of 
statistical  inference  and  decision  theory  within  the  FOT  framework  remains  consistent 
with  that  of  the  stochastic-process:  to  draw  conclusions  with  regard  to  the  infinitely  long 
time-series  from  parameters  measured  from  finite-length  observations  (a  theoretical  subset 
of  the  infinite  series).  Within  this  framework  it  is  possible  to  calculate  bias,  variance  and 
other  measures  of  reliability  and  therefore  pursue  optimizations  as  well  Optimization  of 
an  algorithm  to  estimate  both  time  and  frequency  differences  of  arrival  is  the  core  of  the 
research  presented  here  and  thus  relies  heavily  upon  the  statistical  inference  and  decision 
theory  developed  in  [7],  [14]  and  [15]  for  its  foundation. 

To  that  end,  the  remaining  sections  in  this  Chapter  more  rigorously  develop  the 
cyclostationary  theory  within  the  FOT  model  approach.  In  addition,  to  aid  in  comparison 
to  the  GCC  and  CAF  and  prepare  for  specific  cyclostationary  TDOA  and  FDOA 
algorithms,  digital  implementation  issues  for  cyclostationary  calculations  follow.  This  sets 
the  stage  for  Chapter  IV,  a  survey  of  existing  cyclostationary  TDOA  and  FDOA 
algorithms. 
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B.  CYCLIC  CORRELATION  FUNCTIONS 


1.  Continiious-Tiine  Cyclic  Correlation  Functions 

Traditionally,  a  common  second  order,  time-domain  statistic  used  in  signal 
processing  is  the  autocorrelation  function  given  by  the  quadratic  transformation 

Rxxi'^)-{x{t)x{t-'i:))  (3-13) 

In  the  case  of  a  signal  that  can  be  modeled  as  cyclostationary,  this  transform  will  produce 
spectral  lines  at  non-zero  frequencies  or  such  that 

^  0  (3- 14) 

where 


(3-15) 

This  signal  x{i)  as  noted  above  contains  second  order  periodicity  manifested  in  the 
PSD  of  the  delay  product  given  in  equation  (3-15)  above.  Transforming  (3-15)  into  a 
symmetric  delay  product  and  accommodating  complex  signals  as  well  gives 

y(0  =  (3-16) 


which  forms  the  basis  for  the  fundamental  second  order  cychc  moment  known  as  the  FOT 
continuous-time  (CT)  cyclic  autocorrelation  function: 
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Of  note,  the  CT  cyclic  autocorrelation  function  (drop  the  FOT  as  all  further 
development  remains  within  that  framework)  assumes  the  form  of  the  Fourier  coefficients 
of  the  additive  sine- wave  con^)onents  produced  by  the  periodicity  of  the  delay  product  in 
(3-16).  Also  note  that  (3-17)  is  the  time-average  equivalent  of  (3-7)  from  the  probabilistic 
approach. 

The  CT  cyclic  cross-correlation  function  is  a  natural  extension  of  (3-17)  to  the 
case  of  two  signals.  The  simple  substitution  of  y*(t)  for  /(t)  yields 

(3-18) 

which  serves  as  the  definition  for  the  cyclic  cross-correlation  function. 

Another  interpretation  of  (3-17)  is  as  the  traditional  stationary  autocorrelation 
function  multiplied  by  a  kernel  that  produces  spectral  lines  at  frequencies  where  the 
stationary  autocorrelation  function  contains  additive  sine- wave  components  indicative  of 
its  periodicity.  Consequently,  for  either  interpretation,  at  a  =  0,  the  cyclic  autocorrelation 
function  (drop  the  CT  as  the  argument  apples  equally  to  either  the  continuous  or  discrete 
time  versions  -  see  below)  reduces  to  the  traditional  autocorrelation  function. 

A  final  interpretation  of  the  cyclic  autocorrelation  function  can  be  seen  by  defining 
first 


u{t)  =  x{t)e 


(3-19) 


and 


(3-20) 

so  that  M(r)  and  v(0  are  frequency  shifted  versions  of  x(t).  Now  can  be  written 
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i^(r) = +0 = 


(3-21) 


which  is  the  cross-correlation  of  the  two  versions  of  x(t)  shifted  up  and  down  by 
frequency  ct/2.  In  other  words,  the  cyclic  autocorrelation  function  may  be  viewed  as  the 
correlation  in  the  time  domain  between  two  values  of  x(f)  separated  in  frequency  by  a. 

Consider  as  an  example,  a  generalized  Binary  Phase  Shift  Keyed  (BPSK)  signal, 
not  unlike  the  specific  signal  used  throughout  Chapter  H, 

s(t)  =  Acd{t)  cos[27ifct + 0)  (3-22) 


where  d{t)  is  the  binary  modulating  wave  form  consisting  of  positive  and  negative 
rectangular  pulses  given  by 


eo 

d(f)  df^qif  ~t0~  ) 

n=-oo 


(3-23) 


Approximated  as  a  random  binary  wave,  it  contains  no  first  order  periodicity  and  thus  no 
spectral  lines  in  its  PSD.  Consequently,  the  PSD  of  the  BPSK  signal  j(r)  is  a  scaled  sine 
squared  function  and  will  also  contain  no  spectral  lines,  as  is  evident  its  expression 


sin^cT; 

^cT 


(3-24) 


Multiplying  the  conventional  autocorrelation  function  by  the  cyclic  kernel  the  cyclic 
autocorrelation  ftmetion  follows  as 
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^  J^ct+2jc  /'2^r[or+ 2/cl^o  ^ 


ATu 


(3-25) 


where 

Equation  (3-25)  clearly  shows  additive  sine-wave  components  at  a  =  ±2/.  ±i^and 
a = ±/^[16].  Thus,  the  quadratic  transform  that  is  the  cyclic  autocorrelation  function 
unmasks  hidden  periodicity  in  this  simple  BPSK  signal  and  proves  it  polycyclostationary  in 
the  process. 


2.  Discrete-Time  Cyclic  Correlation  Functions 

As  Gardner  notes  in  [14],  the  concepts  presented  above  apply  transparently  to 
either  the  CT  or  discrete-time  (DT)  case.  Instead  of  the  CT  signal  representations  and  the 
CT-averaging  operation 


Z-^  2Z 


(3-27) 


use  DT  representations  of  the  signals  and  the  DT  counterpart  to  (3-27) 

In  this  manner  it  is  possible  to  produce  the  DT  cyclic  correlation  functions.  For  the 
general  case  of  a  complex  signal,  the  DT  cyclic  autocorrelation  function 
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(3-29) 


rj^Ttat 


which  follows  directly  from  (3-17)  using  (3-28).  Consequently,  the  DT  cyclic  cross- 
correlation  function 

is  the  result  of  a  simple  substitution  of  y*(t)  for  jc*(t). 


C.  SPECTRAL  CORRELATION  DENSITY  FUNCTIONS 


1.  Continuous-Time  Spectral  Correlation  Density  Functions 


The  frequency  domain  equivalent  for  the  CT  cyclic  autocorrelation  function  is  the 
CT  spectral  correlation  density  (SCD)  function.  Gardner  takes  two  approaches  m  [14]  to 
define  the  SCD.  In  the  first,  he  points  out  that  the  SCD  is  simply  an  extension  of  the 
process  by  which  the  power  spectral  density  (PSD)  is  estimated.  This  process  involves 
passing  the  SOI  through  a  bank  of  bandpass  filters,  measuring  the  power  in  each  band  at 
the  output  of  the  filters.  These  filters  constitute  the  traditional  bins  associated  with 
spectrum  analysis.  The  result  of  the  process  is  the  isolation  of  signal  power  in  specific 
frequency  bins.  This  situation  that  may  be  summarized  mathematically  as 


5(/)=lim^ 

b-^oB 


(3-31) 
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where  /if  (r)is  the  impulse  response  of  the  one-sided  bandpass  filter  with  center  frequency 
/and  bandwidth  B. 

The  transition  to  the  SCD  definition  is  a  simple  matter  of  frequency  shifting  the 
SOI  up  and  down  by  a/2,  channelizing  and  correlating.  The  frequency  shifting  operation 
is  the  same  that  occurs  in  the  final  interpretation  of  the  cyclic  correlation  fimction  and 
results  in  uif)  and  v(r)  identical  to  (3-18)  and  (3-19).  Gardner  in  [14]  defines  one  channel 
of  the  SCD  estimation  process  as 

5?(/)=gm-|([*^W®«W][ft^(f)®v(0]*)  (3-32) 

the  time-correlation  of  frequency  shifted,  channelized  segments  of  the  SOI. 

As  a  second  definition,  note  that  is  can  be  shown  that  the  SCD  is  the  Fourier 
Transform  of  the  cyclic  correlation  function  [9].  Consequently  the  SCD  follows  directly 
from  (3-17)  and  has  the  form 


(3-33) 


More  commonly  used  in  cyclostationary  signal  processing,  the  CT  SCD  reduces  to 
the  conventional  PSD  at  a=  0  just  as  the  CT  cyclic  autocorrelation  function  reduced  to 
the  conventional  autocorrelation  function  at  a=  0.  Note  that  (3-33)  represents  the  time- 
average  equivalent  of  (3-8),  the  probabilistic  spectral  correlation  ftuiction. 

In  addition,  the  CT  spectral  cross-correlation  function  follows  directly  from  (3-32) 
and  (3-33)  with  the  appropriate  substitution.  It  has  the  forms 

5«(/)=gm-||ft|(/)®«(r)][ft|(r)®v(f)]*)  (3-34) 


(3-35) 


where  v(f)  now  becomes  a  frequency  shifted  version  of  y(r)  in  (3-34)  and  R^{t)  replaces 
<(r)  in  (3-33)  to  obtain  (3-35). 

Continuing  the  BPSK  exan^)le  from  above,  the  SCD  can  be  found  directly  from 
(3-33).  It  has  the  form 


5?(/)= 4^[e(/+/»  +f)2’(/+/<.-f] 
+d^f-fo  +f )e*(/-/„ 


for  a  =  nlTu  and 


(3-36) 


(/)=  4^[(2(/+/.  +f 

+<^f+fc  +f )e'(/-/„ 


for  a  =  ±Zfg  ±n/Ti,  with  Q(f)  being  the  Fourier  transform  of  the  keying  envelope  g(t) 
[14]. 

The  SCD  is  so  named  as  it  is  in  fact  a  correlation  of  the  SOI  in  the  frequency 
domain.  The  SCD  at  frequency /o  and  cyclic  frequency  is  merely  the  correlation  of  two 
values  of  the  signal  in  the  frequency  domain  separated  in  frequency  by  oj,  and  centered  at 
frequency /o.  Figure  3-1  shows  this  process  concept  graphically. 

The  SCD  is  plotted  on  what  is  called  the  bi-frequency  plane.  The  plane  is  defined 
along  one  axis  as  spectral  frequency /  and  along  the  opposing  axis  as  cyclic  frequency  on 
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Figure  3-1  -  Spectral  Correlation  Function  illustrated,  from  [17] 

Figure  3-2  illustrates  this  point.  The  magnitude  of  the  SCO  corresponds  to  a  height  above 
the  bi-frequency  plane  and  is  often  plotted  in  that  fashion.  Because  or  corresponds  to  the 
separation  distance  between  correlated  values  in  the  frequency  domain  and  the  values  for 

spectral  frequency /naturally  range  from  to  ^^2’  values  for  arrange  from  -/$  to 

u 


2.  Discrete-Time  Spectral  Correlation  Functions 

As  is  the  case  with  the  cyclic  correlation  functions,  the  theory  presented  above  for 
the  CT  SCD  functions  can  be  applied  equally  well  to  the  DT  SCD  functions.  The  DT 
SCD  functions  take  the  form 


s?<j)=  ■  (3-37) 

and 
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(3-38) 


s5(/)= 

'JSS—OO 


Figure  3-2  -  Bi-frequency  plane 


D.  DIGITAL  COMPUTATION  OF  THE  SCD 


1.  Introduction 

The  SCD  computation  can  be  computationally  complex  and  demanding  on  any 
processor.  For  that  reason,  the  development  of  computationally  efficient  algorithms  is 
crucial  to  the  practical  application  of  the  theory.  These  conputational  algorithms  fall  into 
either  frequency-smoothing  or  time-smoothing  categories.  Roberts,  Brown  and  Loomis  in 
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[20]  noted  that  the  frequency-smoothing  algorithms  were  most  efficient  when  estimating 
the  cyclic  spectrum  for  few  values  of  a  or  when  calculating  the  cyclic  spectrum  for  small 
time-frequency  resolution. 

The  time-smoothing  algorithms  prove  more  efficient  for  the  general  case  of 
caicpiaring  the  entire  cyclic  spectrum  of  two  complex  signals.  Specifically,  two  algorithms 
were  developed  in  [20]  with  an  exhaustive  evaluation  of  the  complexity  of  each.  The  Fast 
Fourier  Transform  Accumulation  Method  (FAM)  and  the  Strip  Spectral  Correlation 
Analyzer  (SSCA)  both  simplify  and  thereby  accelerate  the  computation  of  the  SCD. 
Because  cyclostationary  TDOA/FDOA  algorithms  typically  require  some  a  priori 
knowledge  of  the  parameters  of  the  signal  of  interest  and  calculating  the  SCF  of  a  data  set 
directly  from  (3-37)  is  prohibitively  expensive  in  terms  of  computations  required,  digital 
techniques  that  approximate  (3-37)  closely  without  the  computational  burden  are  a 
necessity.  Because  of  this,  FAM  and  SSCA  play  an  in^ortant  role  in  this  work,  providing 
one  method  of  noeasuring  cycle  frequencies  in  the  received  data  sets  without  extensive 
confutation  time.  In  [20],  Roberts,  Brown  and  Loomis  concluded  that  the  SSCA  was  in 
fact  the  most  efficient  of  the  two  methods  for  calculating  the  entire  bi-frequency  plane  and 
thus  it  is  used  extensively  in  the  research  here.  For  that  reason  and  as  an  excellent 
illustration  of  the  compromise  between  performance  and  computational  efficiency,  it  is 
presented  in  more  detail  below. 

2.  Strip  Spectral  Correlation  Analyzer 

Directly  from  its  definition,  (3-37),  the  SCD  can  be  confuted  using  a  traditional  time¬ 
smoothing  approach  as  depicted  in  Figure  3-3.  White  this  approach  provides  the  most 
accurate  estimate  of  the  SCD,  it  proves  prohibitively  time  consuming  and  inappropriate  for 
most  applications  [18]. 

The  SSCA  algorithm,  a  modification  of  this  standard  time-smoothing  method, 
provides  an  excellent  combination  of  accurate  estimation  of  cycle  frequencies  of  interest 
and  computational  efficiency  [18].  By  eliminating  one  of  the  band-pass  filters  in  Figure  3- 
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6  and  using  a  Fourier  transformer  at  the  input  and  output,  SSCA  computes  the  SCD  along 
diagonal  str^s  of  the  bi-frequency  plane  as  depicted  in  Figure  3-4. 

While  the  SSCA  does  manage  to  compute  the  SCD  in  the  most  efficient  manner 
developed  thus  far,  it  makes  some  compromises.  The  output  SNR  suffers  slightly  as  a 
result  of  the  efficiency  that  is  gained  by  the  elimination  of  the  filter  and  the  use  of  the 
Fourier  transform.  In  addition,  the  strips  cause  the  computation  of  the  SCD  to  be  off- 
center  with  respect  to  the  alpha  axis,  thus  distorting  the  gnq)hical  depiction  of  the  SCD. 
However,  the  con^)utational  savings  gained  far  outweigh  the  minor  degradation  in  signal 
power  to  noise  power  at  the  output  of  the  analyzer  and  the  inconvenience  of  the  graphical 
distortion. 

Figure  3-5  illustrates  the  SSCA  architecture.  The  input  filter  consists  of  a  sliding, 
Hamming- windowed,  coarse  FFT  of  lengthiV'.  Each  of  these iV' length  segments  is 
essentially  the  output  of  a  band-pass  filter  with  a  bandwidth  of  1/A/".  This  band-pass 
output  is  then  downconverted  to  baseband  by  a  complex  exponential  multiplication. 

These  products  are  termed  the  complex  demodulates  of  the  input  signal  and  appear  in 


Figure  3-3  -  Time-smoothing  spectral  correlation  analyzer  [19] 
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Figure  3-4  -  SSCA  computation  in  the  bi-frequency  plane,  N'=  8,  from  [19] 


Figure  3-5  -  SSCA  architecture,  from  [19] 


greater  detail  in  Figure  3-6.  As  can  be  seen,  the  process  of  windowing  and  sliding  at  the 
input  filter  essentially  decimates  the  input  sequence.  Brown  and  Loomis  found  that  the 
most  appropriate  value  for  the  decimation  fiictor  or  subsampling  parameter  L  wasAr'/4 
[19]. 


Rgure  3-6  -  Detail  of  the  computation  of  conq)lex  demodulates 

Because  these  complex  demodulates  are  multiplied  by  the  original  sequence  before 
the  final  FFT,  they  must  be  interpolated  to  the  original  sampling  rate.  In  [19]  and  [20] 
Loomis  and  Brown  use  a  simple  hold  operation.  Linear  interpolation  induces  less  high 
frequency  error  and  adds  little  in  computational  complexity.  It  is  used  here  in  place  of  the 
first  order  sample  and  hold  interpolation.  Finally,  the  interpolated  complex  demodulates 
and  the  original  sequence  are  complex  multiplied  and  the  product  Fourier  transformed  by  a 
full  length,  unwindowed  FFT.  The  relationship  of  the  output  to  the  bi-frequency  plane 
appears  in  Figure  3-7.  Each  row  of  output  represents  one  of  theA^' diagonal  strips  in  the 
bi-frequency  plane  as  depicted  in  Figure  3-4.  Each  column  represents  the  values  of  the 
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SCD  along  the  diagonal  lines.  The  resolution  along  lines  of/is  coarse  at  1  /  N'as  a 
byproduct  of  the  computational  savings.  However,  the  resolution  in  a,  at  l/N,  is  much 
higher  which  accommodates  the  characterization  nicely. 


Figure  3-7  -  Block  form  of  SSCA  output 

Given  that  the  majority  of  operations  in  the  SSCA  can  be  in^lemented  as  vector  or 
matrix  operations,  its  implementation  in  any  vector  processing  language  such  as  Matlab® 
is  fairly  sin^)le.  The  psuedocode  for  the  algorithm  jq)pears  in  Figure  3-8.  The  output  of 
the  SSCA  when  implemented  is  a  matrix  withAT'rows  corresponding  to  theiV' diagonal 
strips  of  the  computation  and  columns  corresponding  to  the  SCD  along  the  respective 
str4)s.  Referring  again  to  Figures  3-4  and  3-6,  these  str^s  are  not  confuted  along  lines  of 
/and  or  and  must  therefore  be  transposed  along  those  lines  before  plotting. 

The  transformation  involves  applying  the  relationship  between  the  diagonal  strip 
subscript  k,  the  frequency/ and  the  cyclic  frequency  a  Those  relationships  from  [20]  are 
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(3-39) 


—AT'  N' 

-^<k<— — 1 
2  2 


«  =  2/^-2/ 


(3-40) 


While  cyclostationary  signal  processing  techniques  provide  powerful  tools  for  use 
in  poor  signal  environments,  those  tools  are  used  at  the  expense  of  computational 
coiiq)lexity.  The  SSCA  strives  to  reduce  that  complexity  without  serious  degradation  in 
the  SNR  of  the  output.  In  the  next  chapter,  specific  applications  of  cyclostationary  signal- 
selectivity  to  the  estimation  of  TDOA  are  surveyed.  In  this  section,  it  will  become  even 
more  clear  the  value  that  computationally  efficient  digital  implementations  have  in  the 
practical  en^loyment  of  some  of  the  cyclostationary  techniques. 


/*  Compute  Complex  Demodulates  of  x  */ 

Do  p:ss  0  to  P-1 

Compute  =  FFT[a(r)x(pI+r)] 

Dok:=-Ny2toN'/2-l 

Compute  =  xT(pL,fj^)&xp{-j2^Lk/N‘ } 

end 
end 

/*  Interpolate  Xj-fpL,^)  */  ^ 

Compute Xj(pL, Xj.(nyf^) 

/*  Compute  and  Smooth  Product  Waveforms 
Doife:=  -N'/2toN'/2-\ 

Compute  SiM>fk)  =  y\n)Xr(n,f^) 

Compute 

end 


Figure  3-8  -  Psuedocode  for  SSCA  computation,  from  [19] 
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IV.  CYCLOSTATIONARY  TDOA  AND  FDOA  ESTIMATION 


A.  CYCLOSTATIONARY  TIME  DIFFERENCE  OF  ARRIVAL 


1.  Cyclostationaiy  TDOA  Signal  Model 

The  signal  model  used  to  develop  cyclostationary  TDOA  algorithms  is  Model  B 
and  identical  to  that  presented  in  Chapter  H.  Recall  the  signals  received  at  two  spatially 
separated  observers  can  be  given  by  (2-1)  and  (2-2)  and  are  repeated  below  for  clarity 

x(t)  =  s(t) + nj  (0  (4-1) 

and 

y(t)  =  Ar-s(t-D)+n2(t)  (4-2) 

The  difference  between  the  two  models  lies  in  the  temporal  and  spectral  relationships 
between  s(t),  ni(t)  and  n2(t).  In  Chapter  H,  ni(t)  and  nzO)  were  assumed  to  have  no 
con^onents  that  temporally  and  spectrally  overlap  the  SOI,  ^(t).  In  the  cyclostationary 
TDOA  signal  model,  ni(t)  and  nzit)  represent  all  SNOI(s)  and  noise  present  at  the 
respective  observers.  Tl»y  may  or  may  not  contain  co-channel  interferers  in  reality, 
however  for  purposes  of  this  work,  they  are  assumed  to  contain  interference  that  in  fact 
q)ectrally  overlaps  the  SOL 

Given  this  model,  the  cyclic  auto  and  cross-correlation  functions  in  general  can  be 
written  as 


/l?(T)=Jif(T)+Jl^(T)  (4-3) 

(4-t) 
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(4-5) 


just  as  (2-3)  -  (2-5)  represented  the  conventional  auto  and  cross-correlation  functions.  In 
addition,  as  with  (2-6)  -  (2-8),  the  spectral  correlation  functions  are 

(4-6) 

^(J)  =  \ArfsS’(f)e->'‘’^+S^(J)  (4-7) 

=  (4-8) 


The  cross- SCD(s)  for  the  SOI  and  the  noise  are  omitted  because  they  are 
identically  zero  with  the  assumption  that  the  SOI  and  the  noise  are  independent.  This 
assumption  is  not  made  for  any  SNOI  contained  in  the  noise  terms.  They  are  assumed  to 
be  correlated  and  thus  the  cross-SCD  for  the  noise  terms  at  each  receiver  appears  in  (4-8). 

Though  (4-6)  -  (4-8)  contain  the  TDOA,  D,  this  parameter  is  now  ma^d  by  the 
spectrally  overlapping  conq)onents  of  ni(t)  and  n2(r).  Thus,  any  atten^t  to  estimate  D 
with  traditional  methods  that  compute  the  above  equations  at  or  =  0,  like  those  illustrated 
in  Chapter  n,  would  result  in  a  corrupted  value.  However,  if  j(f)  contains  some  cycle 
frequency,  Oo,  for  which  the  cyclic  cross-correlation  between  ni(f)  and  nzit)  is  zero,  then 
cyclostationary  techniques  can  determine  the  TDOA,  D  of  the  SOI,  5(r)  even  in  the 
presence  of  strong  contributions  from  ni(t)  and  nz(,t)  that  would  otherwise  corrupt  the 
estimate.  Essentially  this  eliminates  the  SNOI(s),  resulting  in  (4-6),  (4-7)  and  (4-8)  being 
rewritten  for  the  use  of  oj)  as 


5?»(/)=Sf”(/) 

(4-9) 

(4-10) 
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Sy^  (/) = ArSf‘> 


(4-11) 


Thus  a  reliable  estimate  of  D  is  possible  even  in  highly  corrupt  environn^nts 
provided  a  unique  Ok  exists  for  the  SOL 

2.  Current  Cyclostationaiy  TDOA  Determination  Algorithms 

In  [16],  Gardner  and  Chen  outlined  no  less  than  seven  different  cyclostationary 
TDOA  algorithms.  The  algorithms  range  in  complexity  from  simple  cyclic  modifications 
to  the  GCC  to  more  complicated  minimum  least  squares  approaches,  each  with  its 
strengths  and  weaknesses.  The  three  algorithms  of  this  group  that  most  thoroughly 
represent  the  range  of  techniques  developed  in  [16]  are  described  below. 

The  Cyclic  Cross-Correlation  (CCC)  method  is  the  result  of  the  simple  substitution 
of  the  cyclic  cross-correlation  function  for  the  limit  auto-correlation  function  in  the  CC 
TDOA  method  from  above,  yielding  the  following  relation  from  (49)  in  [16] 

a„(r)^i^,(T)|=|Aj/^(T-D)|  (4-12) 

Note  that  the  inclusion  of  the  cyclic  cross-correlation  function  in  the  solution  injects  the 
inherent  signal  selectivity  of  cyclostationary  signal  processing  into  the  CC-like  TDOA 
determination,  significantly  mitigating  the  nagging  problem  of  co-channel,  co-temporal 
interference  outlined  with  the  CC  and  CAF  above.  Gardner  and  Chen  found  that  although 
the  peak  of  («)|  does  not  occur  at  u  =  0  for  all  signals,  since  it  is  an  even  function  of  u, 

in  those  situations  lacking  a  single  peak  at  u  =  0,  usually  a  pair  of  peaks  straddle  u  =  0 
thereby  identifying  the  TDOA  unambiguously  [16]. 

Gardner  and  Chen  further  point  out  that  the  performance  of  the  CCC  can  be 
significantly  improved  using  cyclic  versions  of  the  same  methods  employed  in  the 
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construction  of  the  GCC.  The  direct  cyclic  counterpart  of  the  GCC  in  (2-15)  can  be 
written  as  the  estimate  bjT)of  the  ideal  measurement 


(4-13) 


which  is  (56a)  in  [16]  and  always  produces  a  peak  at  the  TDOA  value.  Because  it  is 
essentially  the  ratio  of  cyclic  cross  spectrum  to  the  cyclic  auto-spectrum  at  a  selected  cycle 
frequency,  this  method  is  known  as  the  Spectral  Correlation  Ratio  (SPECCORR), 

A  significant  difference  between  this  cyclic  version  of  the  GCC  and  the  traditional 
GCC  is  the  phase  content  of  the  weighting  functions  used  in  SPECCORR.  The  weighting 
functions  seen  in  Table  2-1  above  merely  emphasize  or  de-emphasize  certain  spectral 
bands  essentially  operating  effectively  only  at  a  =  0..  The  weighting  functions  in 
SPECCORR  also  contribute  to  the  band  selection  along  that  a=0  line  in  the  cyclic  plane. 
However,  the  phase  content  contained  in  the  complex  values  of  these  cyclic  weighting 
functions  also  serve  to  pre-filter  further  in  the  cyclic  plane  to  eliminate  interference  and 
emphasize  features  at  cycle  frequencies  other  than  cr=  0.  Herein  lies  the  performance  gain 
The  third  estimator  in  [16],  Spectral  Coherence  Alignment  (SPECCOA)  was 
developed  and  presented  as  the  solution  to  an  ad  hoc  minimum  least  squares  (MLS) 
approach.  Because  it  has  proven  an  excellent  performer  in  a  broad  range  of  environments, 
more  emphasis  will  be  placed  on  its  development  here. 

Using  equations  (4-3)  -  (4-5),  and  the  assurtq)tion  that  s(t)  contains  a  unique  Ok,  it 
can  be  seen  that 


R“{u)=CR“(,u-D)e-i’^  (4-14) 

The  cross-correlation  term  for  ni(jt)  and  mit)  is  eliminated  by  the  use  of  a  cycle  frequency 
unique  to  the  SOI.  An  estimate  of  the  TDOA,  D  that  minimizes  the  magnitude  of  the  sum 
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of  the  square  errors  between  the  measured  value  of  the  left  side  of  (4-14),  and  the 

measurement  of  the  right  side  of  (4-14),  (with  D  substituted)  yields  the  minimum- 
least-squares-optimized  value  for  the  TDOA. 

Mathematically,  from  [16], 


£)=argmin 

A 

c,r 


(4-15) 


where  Ca(T)  is  the  estimate  of  the  MLS  fiinction.  It  was  further  shown  in  [16],  that  the 
solution  of  the  MLS  problem  is  of  the  form 


where 


Z)=argmax 


(4-16) 


(4-17) 


(4-18) 


where  (4-18)  derives  fi’om  (4-17)  through  Parseval’s  relation.  Gardner  and  Chen  go  on  to 
further  prove  that  (4-18)  does  indeed  peak  at  t=  D. 

SPECCOA  derives  its  name  from  the  fact  that  the  peak  in  (4-18)  occurs  by 
maximizing  the  correlation  in/ of  the  two  spectral  correlation  functions  through  the 
alignment  of  the  phases  of  the  respective  functions. 

In  [4],  Spooner  and  Gardner  derive  two  additional  cyclostationary  TDOA 
algorithms,  both  considerably  more  computationally  complex  than  the  three  outlined 
above.  The  first  is  a  more  extensive  generalized  version  of  SPECCOA  termed  G- 
SPECCOA  and  the  second.  Maximum  Likelihood  (ML),  is  the  result  of  a  maximum 
likelihood  development. 
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The  generalization  of  G-SPECCOA  adds  seven  additional  terms  to  the  estimator 
taking  into  account  multiple  cycle  frequencies  and  both  conjugate  and  non-conjugate 
cyclostationarity.  It  includes  the  products  of  the  cross-spectral  correlation  and  the  auto 
spectral  correlation  which  are  referred  to  as  the  auto-cross  terms.  In  addition,  G- 
SPECCOA  includes  products  of  cross  spectral  correlations  and  cross  spectral  correlations 
(cross-cross  terms)  and  products  of  auto  spectral  correlations  and  auto  spectral 
correlations  (auto-auto  terms).  The  result  is  an  eight  term  algorithm  that  has  the  form 


argmax 

d 


Re 


a 


— oo 

— OO 

— oo 


(4-19) 


Finally,  Spooner  and  Gardner’s  ML  algorithm  is  derived  for  Model  B  modified  to 
include  complex  constants  received  signals.  Through  a  complex  series  of  steps  and 
simplifications,  some  of  which  are  utilized  in  the  ML  derivation  in  Chapter  V,  Spooner  and 
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Gardner  derive  an  asynq)totically  ML  estimator  for  TDOA  that  also  has  eight  terms  and  a 
computational  complexity  similar  to  that  of  G-SPECCOA. 

Like  G-SPECCOA,  the  algorithm  includes  auto-cross  terms,  cross-cross  terms  and 
auto-auto  terms.  In  its  ideal  ML  form,  the  algorithm  is  expressed  as  a  function  of  the 
receiver  data  records  and  the  SOI  waveform.  As  complete  knowledge  of  the  SOI  is 
usually  not  available,  this  algorithm  is  considered  an  ideal  ML  TDOA  estimator. 
Substituting  for  the  terms  that  contain  the  SOI  with  approximations  that  contain  only  the 
data  records,  converts  this  ideal  algorithm  into  one  that  is  implementable.  An  infinite 
number  of  substitutions  are  possible  so  Spooner  and  Gardner  propose  several 
approximations  in  [4]  to  create  practical  estimators.  The  original  ML  algorithm,  equation 
(101)  in  [4]  appears  below  as  (4-20). 
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(4-20) 
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V. 


NEW  MAXIMUM-LIKELfflOOD  ESTIMATORS 


A.  BACKGROUND 

1.  Signal  Modd 

The  ultimate  goal  of  this  work  is  to  find  the  ML  estimators  for  the  TDOA  and 
FDOA  of  a  weak  cyclostationary  signal  in  white  Gaussian  noise.  To  that  end,  an 
appropriate  signal  model  must  be  used  both  to  accurately  reflect  the  signal  environment  of 
interest  and  to  make  the  derivation  mathematically  tractable.  Signal  Model  C  is  chosen  to 
most  accurately  depict  the  case  of  two  receivers  collecting  a  single  SOI  in  white  Gaussian 
noise  while  maintaining  some  degree  of  mathematical  tractability.  One  change  made  to  the 
model  is  the  elimination  of  the  arbitrary  gain  and  phase  mismatch  between  the  receivers. 
This  complex  constant  scaling  the  output  of  receiver  #2  represents  a  mismatch  that  is  an 
unnecessary  distraction  in  this  derivation  for  several  reasons.  First,  gain  mismatch  is  easily 
accommodated  with  some  form  of  normalization  if  such  a  procedure  is  even  necessary. 
Second,  modem  collections  systems  enjoy  the  benefits  of  digital  technology,  GPS  and  fine 
frequency  standards  and  thus  take  advantage  of  excellent  coherence  in  both  timing  and 
frequency.  Finally,  the  inclusion  of  such  a  complex  constant  in  a  maximum  likelihood 
derivation  would  involve  a  maximum  likelihood  estimate  of  its  own,  taking  away  from  the 
original  intent  of  the  research.  Thus,  to  concentrate  on  a  joint  ML  estimate  for  TDOA  and 
FDOA,  the  gain  and  phase  mismatch  can  be  dropped  here. 

The  new  model  for  the  derivation  of  the  MLE  is  a  sin^)lified  Doppler  signal  model 
similar  to  Model  C  but  expressed  in  discrete  time  and  lacking  the  complex  constant  in  the 
second  receiver  equation: 


xin)—s{n)+mx{n) 


(5-1) 
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+my{n) 


(5-2) 


y{n)=s{n-T>)e 


+j2frf^n 


for  »  =  1, 2, A^and  a]l  but  n,  D  and  fd  are  complex- valued. 

2.  Matiiematical  Model 

To  find  ML  estimates  for  D  and  fi,  it  is  necessary  to  maximizf^.  a  likelihood 
function  conditioned  on  only  those  parameters.  In  order  to  do  so,  a  likelihood  function  for 
an  equivalent  complex- valued  model  conditioned  on  all  the  unknown  parameters, 
including  but  not  hmited  to  D  and  ^  can  be  formed  [4].  In  due  course,  the  conditioning 
on  the  unknown  parameters  not  of  interest  are  removed  to  leave  the  likelihood  function 
needed  for  the  estimates  of  D  and  fj.  The  equivalent  complex- valued  model  is 

z  =  y+w  (5-3) 


where 


'x(n) 

s(n) 

y(n) 

j  v  = 

and 


w  = 


mx(n) 

ntyin) 


(5-4) 


and  the  dimensions  of  each  vector  set  at  2Ar  x  1. 

For  the  conditioning  on  the  unknown  parameters,  let  b  and  0  denote  the  random 
and  nonrandom  conditioning  parameters  respectively.  To  find  the  likelihood  function,  the 
probability  density  function  (pdf)  of  z  is  needed  specifically  conditioned  on  b  and  0  or 
equivalently,  conditioned  on  (p  with 


^=[b  -  (5-5) 

i 

I 
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where 


e=\p  fj]  (5-6) 

Thus,  the  equivalent  complex  valued  mathematical  model  and  the  vector  of 
conditioning  parameters  contain  all  the  information  about  the  signal  model  needed  to 
develop  a  likelihood  function  that  will  produce  maximum  likelihood  estimates  for  the  two 
parameters  of  interest  D  and  fd. 

B.  DERIVATION 

Of  key  importance  in  the  formation  of  this  likelihood  function  is  the  statistics  of  the 
noise.  First,  the  only  knowledge  of  the  system’s  statistics  is  the  definition  of  the  noise  at 
each  receiver  as  zero-mean,  additive  white  Gaussian  noise.  Second,  the  noise  terms  are 
assumed  to  be  uncorrelated  and  thus  to  contain  not  SNOI.  Thus  the  noise  pdf, 
conditioned  on  the  signal  parameters  completely  describes  the  system.  Expressing  the 
noise,  /nx(n)  and  %(«),  contained  in  the  vector  w,  in  terms  of  its  signal  and  receiver 
conponents  produces  a  pdf  with  receiver  terms,  x{ri)  and  y(/i)  contained  in  the  vector  z, 
and  signal  terms,  sin)  and  its  time-delayed  and  shifted  version  contained  in  the  vector  v, 
conditioned  on  all  the  unknown  random  and  non-random  parameters.  It  takes  the  form  of 
a  discrete  jointly  Gaussian  pdf 

|exp[-(z-v)*^(R,_,)‘‘(*-v)j  (5-7) 

where  againw  =  z— v  [21]. 

Now  because  the  noise  is  zero-mean  Gaussian,  some  simplifications  can  be  made. 
Specifically,  substituting 


(5-8) 
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and 


(Rz-v)  -J^^NxlN  (5-9) 

produces  the  followmg  changes  to  the  likelihood  function, 

(5-10) 

which  in  turn  simplifies  to 

Finally,  when  expanded  within  the  exponential,  the  expression  becomes 

Pw(w|<P)=^^^^^^;yeXp  -^(zT2-t’Tz-z‘’'v  +  t’Tv)  (5-12) 

Though  highly  relevant  to  an  hardware  and  system  implementation  problem  where 
thresholds  and  amplifier  gains  are  major  issues,  the  constant  of  proportionality  does  not 
depend  upon  the  variables  of  concern  here  and  is  merely  a  DC  offset  in  any  effort  to 
maximize  an  equation  derived  from  the  remainder  of  the  expression  above.  Thus  it  is 
dropped  for  convenience.  Hence,  the  likelihood  function  now  becomes, 

^  ^  *T  ^  ^  *T  \  (  *T  ^ 

P.Hp)o=exp  exp  exp  ^  exp  (5-13) 
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The  weak  signal  assumption  allows  the  expansion  of  the  exponential  expressions 
into  Taylor  series  as  follows 
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(5-14) 


(5-15) 


(5-16) 


(5-17) 


and  substituting  the  four  expressions  into  the  pdf  gives  us  a  straight-forward  but  tedious 
problem. 
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(5-18) 


Expanding  the  quadratic  and  retaining  only  fourth  or  lower  order  terms,  the 
expression  becomes  a  sum  of  fourteen  products 
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(5-19) 


which  are  numbered  Pi  to  Pm  (for  product  1,  product  2, ...)  from  left  to  right,  top  to 
bottom,  discounting  the  leading  constant. 

Without  loss  of  generality,  if  the  sample  rate  is  greater  than  the  Nyquist  rate  to 
ensure  perfect  reconstruction  of  the  analog  signal,  the  transformation  of  the  fourteen 
expressions  above  into  continuous-time  integrals  is  possible.  The  terms  now  are 
numbered  CPi  to  CPm,  for  continuous-time  products  1  to  14.  At  the  same  time,  the 
substitution  of  the  original  system  and  signal  terms  back  into  the  expressions  gives  us  an 
initial  insight  into  the  content  of  the  terms.  To  avoid  the  additional  clutter  associated  with 
the  addition  of  the  time-delay  variable  D  and  frequency  shift  complex  exponential,  the 
variable  Si(t)  represents  the  time-delayed  and  frequency  shifted  version  of  pure  signal  that 
is  measured  at  receiver  2.  The  transformed  terms  appear  below: 

^  fTT  TT 

CP.:  JJ/WxWdr+JJ/WKf)* 


(5-20) 
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Next,  the  weak  signal  assumption  eliminates  term  CP2  because  contribution  to  the 
MLE  of  TDOA  and  FDOA  for  this  term  will  be  very  small  given  that  the  signal  content  in 
term  CP2  is  small  compared  to  the  noise  content.  The  same  is  true  for  the  contribution  of 
the  fourth-order  signal  term  CP14.  It  too  can  be  dropped  because  the  weak  signal 
assumption  causes  fourth-order  terms  to  become  very  small  [4]. 

Finally,  writing  terms  CPi  and  CP13  as 


reveal  they  are  essentially  power  terms.  If  r  is  sufficiently  large,  these  terms  converge  to 
scaled  average  powers  of  the  received  waveforms,  the  SOI  and  the  delayed  and  Doppler- 
shifted  SOI.  These  average  powers  do  not  depend  on  the  TDOA,  D  or  the  FDOA,  fd. 
Consequently,  they  are  dropped  as  contributing  little  to  the  MLE  [4]. 

Following  that  simplification,  the  next  step  should  remove  the  conditioning  on  the 
random  parameter  vector  b.  For  nearly  all  digital  communications  signals,  the  only 
contribution  to  this  random  conditioning  vector  comes  from  the  bit  or  symbol  sequence  of 
the  message  [4].  To  that  end,  if  the  expected  value  of  the  signal  s(t)  is  zero,  as  is  the  case 
in  almost  every  man-made  digital  signal,  carrying  out  the  expectation  over  b  causes  any 
term  with  an  odd  power  of  signal  con^onents  to  go  identically  to  zero.  This  eliminates 
terms  CP3,  CP5,  CPe,  CP9,  CPio  and  CPn. 

This  leaves  terms  CP4,  CP?,  CPs  and  CP12  which  will  collapse  using  the  complex 
conjugate  relationships,  CP12  =  CP/  and  CP?  =  CPs*,  into  two  new  terms  T?  and  T2  with 
the  following  relationship  to  the  original  14  CP  terms: 


Ti  =  CP4  +  CP12  =  2  Re{CP4} 

(5-36) 

T2  =CP7  H-CPg  =  2Re{CPg} 

(5-37) 

with  the  constant  out  front  dropped  as  merely  a  bias  term.  The  new  terms  are  now 
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Note  that  the  first  integral  in  Ti  does  not  contain  either  parameter  of  interest.  The 
same  is  true  for  the  first  integral  in  T2.  Therefore,  these  two  integrals  are  removed  from 
the  estimator.  Now,  carrying  out  the  expectation  with  the  remaining  terms,  enables  the 
expression  of  the  sum  of  products  in  terms  of  the  cyclic  correlation  functions.  The 
expectation  operator  allows  the  treatment  of  the  signal  terms  as  random  variables  and  the 
receiver  terms  as  non-random  constants,  forming  the  probabilistic  auto-  and  cross¬ 
correlation  functions  according  to 

^pq[h^h)=^p(h)^*[h)}  (5-40) 


where  pi-)  and  ^(0  each  represent  one  of  the  signals  {j(  ),  j*(  ),  j^O),  j^(*)}- 

Following  that,  the  remaining  pieces  comprise  the  finite-time-averaged  auto-  and 
cross-correlation  functions  for  various  combinations  of  the  received  waveforms,  for 
example. 
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with  u{t)  defined  as  the  unit  step  function 


otherwise 


(5-43) 


(5-44) 


Finally,  the  cyclic  correlation  functions  follow  using  the  following  relationship 

below 


=  (5-45) 

a 

where  p{t)  and  q(t)  again  represent  one  of  the  signals,  {^O),  /(•)» 

Thus,  the  likelihood  function  may  now  be  expressed  as  products  of  the  cyclic  auto- 
and  cross-correlation  functions  of  the  receiver  and  signal  variables.  The  expression  now 
has  the  form 
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Next,  a  final  substitution  of  sit)  for  Siit)  expresses  the  estimator  in  terms  of  the 
received  waveforms  x(f)  and  yit)  and  the  limit  cyclic  correlation  functions  of  sit).  This  is 
the  Ideal  Two-Dimensional  Maximum  Likelihood  Estimator  (I2DMLE). 
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where  the  expression  to  be  maximized  simultaneously  over  D  and  fd  to  obtain 
cyclostationary  ML  estimates  for  TDOA  and  FDOA. 

Note  that  I2DMLE  is  in  fact  ideal  because  it  requires  knowledge  of  the  ideal  cyclic 
auto-correlation  function  of  the  SOI  a  priori  In  practice  this  cyclic  correlation  function  is 
not  available  thus  requiring  the  ideal  caveat.  Nonetheless,  this  is  the  first  of  three 
cyclostationary  TDOA/FDOA  estimators  implemented  and  tested  here.  The  performance 
of  this  ideal  algorithm  is  instructive  in  assessing  the  loss  in  performance  due  to  the 
approximations  made  in  order  to  realize  the  next  two  algorithms  outlined  below.  In 
addition,  should  an  acceptable  extension  of  the  Cramer-Rao  Lower  Bound  be  found  for 
the  cyclostationary  TDOA/FDOA  case,  the  extent  to  which  the  ideal  algorithm  approaches 
this  bound  as  the  collect  length  increases  will  be  a  valuable  investigation  as  well 

Finally,  in  order  to  more  easily  implement  I2DMLE,  the  expression  can  be  written 
in  terms  of  the  spectral  correlation  functions  by  using  Parseval’s  relation.  I2DMLE  in  the 
frequency  domain,  still  with  maximization  over  D  and  fd,  assumes  the  form 
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C.  ALGORITHMS 


1.  Approximated  2-Dimeiisioiial  Maximum  Likelihood  Algorithm  #1 

To  obtain  the  final  two  estimators,  both  approximations  to  the  I2DM1J5,  we  must 
substitute  for  the  ideal  cyclic  auto-correlation  function  of  the  SOI  with  the  cyclic 
correlation  functions  of  the  data  record  x(t)  and  y(t).  This  is  possible  and  represents  a 
plausible  approximation  because  by  modeling  the  received  signals  as  cyclostationary  and 
choosing  SOl-unique  cycle  frequencies,  in  the  limit  as  T  -»  «>,  the  SOI  spectral  correlation 
function  can  be  measured  as  accurately  as  desired.  Note  as  an  example,  that  the  following 
is  true 


h;mR“(r)  =  R“(T),  a^O 


(5-49) 


and 


BmR|,(T)  =  R£(r)  (5-50) 

Consequently,  the  simplest  practical  approximation  to  the  ideal  estimator  is 
obtained  by  substituting  x(t)  directly  for  s(t)  in  the  above  expression  to  obtain  the 
Approximated  Two-dimensional  Maximum  Likelihood  Estimator  #1  (A2DMLE#1)  which 
is  maximized  over  D  and  ^  and  takes  the  form 
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(5-51) 


Note  that  A2DMLE#1  contains  both  conjugate  and  non-conjugate  terms  with 
cross-cross  terms  and  auto-cross  terms,  resembling  some  portions  of  G-SPECCOA. 
Thus,  it  is  expected  that  the  resulting  combination  possess  some  of  the  strengths  and 
weaknesses  of  this  previously  derived  algorithm  in  some  form. 

In  the  frequency  domain,  A2DMLE#1  looks  identical  to  I2DMLE  withx(0 
substituted  for  s{t)  and  still  maximized  over  D  and  fd. 
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(5-52) 


2.  Approidiiiated  2-Dimeiisioiial  Maximum  Likelihood  Algorithm  #2 

The  second  approximation  to  the  I2DMLE  is  a  logical  extension  of  A2DMLE#1 
above  with  three  concurrent  goals:  (1)  keep  the  strongest  terms  to  maintain  as  much  as 
possible  the  performance  of  A2DMLE#1,  (2)  maintain  both  the  conjugate  and  non¬ 
conjugate  cyclostationary  exploitation  and  (3)  minimize  the  computational  burden  by 
limiting  the  number  of  terms  as  much  as  possible.  Observing  first  that  the  final  summation 
depends  little  on  the  frequency  parameter  of  interest  fd  it  can  be  eliminated  from  the 
second  approximation.  Also  noting  that  the  first  and  second  summations,  as  well  as  the 
fourth  and  fifth  summations,  differ  only  in  their  first  cross-correlation  terms,  only  one  of 
each  of  these  pairs  of  terms  is  retained.  Finally,  to  obtain  another  auto-cross  term  like 
SPECCOA,  the  third  sum  is  converted  into  such  a  term  by  substituting  for  R®  (T)with  the 
following  cyclic  cross-correlation 


KiT^) = (r- 


(5-53) 
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producing  an  auto-cross  term  which  is  theoretically  a  strong  contributor  to  the  MLE . 

The  result  is  Approximated  Two-dimensional  Maximum  Likelihood  Estimator  #2 
(A2DMLE#2),  a  three  term  algorithm  whose  performance  should  not  deteriorate 
considerably  from  that  of  A2DMLE#1,  still  exploits  both  conjugate  and  non-conjugate 
cyclostationarity  and  is  much  more  easily  implemented  with  less  computational  burden.  It 
retains  the  need  for  maximization  over  D  and  fd  and  has  the  form 
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(5-54) 


In  the  frequency  domain,  A2DMLE#2  is  similar  to  A2DMLE#1  with  the  exception 
of  the  second  term  which  has  the  cyclic  cross-correlation  substituted  in  place  of  the  cyclic 
auto-correlation  function.  The  entire  expression  in  the  frequency  domain  has  the  form 
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Now  that  three  different  original  ML  algorithms  are  derived,  a  method  with  which 
to  measure  and  compare  their  performance  is  needed.  In  the  next  chapter,  the  details  of 
the  Monte  Carlo  testing  provide  a  measure  of  the  Root  Mean  Square  Error  in  TDOA  and 
FDOA  estimation  for  each  of  the  three  algorithms.  This  performance  measure  is 
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calculated  for  a  variety  of  signal  environments  and  then  compared  to  the  performance  of 
traditional  techniques  and  previously  derived  cyclostationary  TDOA  algorithms. 
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VL  MONTE  CARLO  PERFORMANCE  TESTING 


A.  OBJECTIVES 

The  Monte  Carlo  Testing  (MCT)  satisfies  four  objectives,  two  quantitative  and 
two  qualitative.  The  quantitative  objectives  for  the  MCT  are:  (1)  provide  a  numerical 
measure  of  the  performance  of  the  three  newly  derived  algorithms,  I2DMLE,  A2DMLE#1 
and  A2DMLE#2  and  (2)  provide  a  quantitative  means  of  direct  comparison  to  the 
performance  of  the  CAP  under  identical  circumstances.  The  qualitative  objectives  are:  (1) 
provide  a  qualitative  means  of  comparison  to  the  performance  of  previously  derived 
cyclostationary  TDOA  algorithms;  and  (2)  provide  a  means  of  comparison  to  the 
stationary  CRB. 

Satisfying  the  first  quantitative  objective  by  empirically  measuring  the  performance 
of  the  three  algorithms  allows  specific  conclusions  to  be  drawn  regarding  the  loss  in 
performance  due  to  each  of  the  two  approximations  to  the  ideal  algorithm.  This  in  turn 
can  lead  to  other  more  informed  approximations  in  future  work.  The  second  quantitative 
objective  allows  the  direct  numerical  con^arison  of  three  cyclostationary  TDOA-FDOA 
estimators  to  the  traditional,  well-known  CAP  under  the  exact  same  conditions.  This 
reveals  to  some  degree  the  extent  of  the  performance  gain  achieved  by  these 
cyclostationary  techniques  and  also  to  some  degree,  under  which  circumstances 
cyclostationary  techniques  are  superior  or  inferior  to  the  traditional  TDOA-PDOA 
estimation  process. 

The  first  qualitative  objective  is  satisfied  by  intelligent  conq)arison  of  the  MCT 
results  with  those  obtained  in  previously  published  work  under  similar  but  not  identical 
testing  scenarios.  The  relationship  of  the  performance  of  these  new  estimators  to  the 
performance  of  those  estimators  that  constitute  the  heritage  for  this  work  is  m^ortant  in 
understanding  if  in  fact  the  TDOA  estimates  found  here  compare  well  with  that  previous 
work.  Additionally,  whether  those  estimates  for  TDOA  were  obtained  in  conjunction  with 
comparable  performance  in  PDOA  estimation  is  an  important  consideration. 
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The  second  qualitative  objective  comes  in  the  comparison  of  the  stationary  CRB 
on  the  performance  of  the  CAP  to  the  results  of  the  MCT.  Since  the  CRB  represents  the 
theoretical  lower  limit  on  the  variance  of  the  estimates  obtained  from  the  CAP,  it  is  in  a 
sense,  the  very  best  that  a  traditional  TDOA-PDOA  estimator  can  do.  Though  the  CRB 
used  here  was  derived  using  different  assumptions  than  those  that  form  the  basis  for  the 
cyclostationary  TDOA-PDOA  estimators  developed  here,  comparing  the  MCT  results  to 
its  ultimate  lower  limit  provides  another  piece  of  evidence  in  support  of  the  added 
contribution  that  cyclostationary  techniques  make  over  and  beyond  even  the  ideal 
stationary  estimator. 

B.  TEST  ENVIRONMENTS 


1.  Introduction 

In  each  case,  the  SOI  is  the  BPSK  signal  introduced  in  Chapter  n  and  detailed  in 
the  Appendix.  Because  this  is  the  only  modulation  scheme  investigated  here  it  hardly 
represents  an  exhaustive  examination  of  the  performance  of  these  algorithms.  But,  BPSK 
is  a  rich  modulation  scheme  from  a  cyclostationary  perspective,  exhibiting  both  conjugate 
and  non-conjugate  features  as  well  as  features  associated  with  the  carrier  frequency,  the 
bit  rate  and  several  combinations  of  the  two.  Thus,  it  is  considered  an  excellent  subject 
for  initial  investigation  with  other  modulation  schemes  saved  for  follow-on  work. 

Regardless  of  modulation  scheme,  three  simple  parameters  dominate  the 
performance  of  TDOA  and  PDOA  estimators:  the  type  of  interference,  the  relative 
strength  of  that  interference  and  the  amount  of  data  available  to  process.  In  a  benign 
situation,  slightly  overlapping  or  no  interference  with  high  SINR  and  long  collect  tengths 
represent  the  best  situation  an  estimator  could  encounter.  On  the  contrary,  completely  co¬ 
channel  and  co-temporal,  partially  co-spatial  interference  with  low  SINR  and  short  collect 
lengths  represent  the  more  realistic  and  decidedly  poor  situation.  Therefore,  in  order  to 
properly  exercise  these  algorithms  in  pursuit  of  the  first  two  objectives  and  allow  for  easy 
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comparison  to  previous  work,  two  of  the  most  stressing  of  environments  are  used: 
wideband  interference  and  narrowband  interference.  These  two  environments  represent 
similar  scenarios  used  in  [4]  and  [22]  thus  facilitating  the  comparison  of  the  results. 

2.  Wideband  Interference  Environment 

Environment  I  tests  the  algorithms  in  a  wideband  interference  environment.  In  this 
case,  the  bandwidth  of  the  SNOI  is  nearly  twice  that  of  the  SOI,  ensuring  complete 
spectral  and  temporal  overlap  as  shown  in  Figure  6-1.  Note  that  the  presence  of  a  second 
signal  is  not  obvious  in  this  plot  of  the  spectral  density  of  the  received  signal  for  receiver 
#1.  In  fact,  the  BPSK  spectrum  that  appears  in  the  plot  is  primarily  due  to  the  SNOI,  not 
the  SOI,  with  the  only  evidence  of  the  SOI  being  the  small  amount  of  additional  energy 
that  causes  the  slight  spreading  toward  the  origin  of  the  SNOI  main  lobe. 


800i  1  i  i'll'  I  '  I  I 


^-4-3  -2  -1  0  1  2  3  4.5 

frequency  (Hertz)  ^ 


Figure  6-1  -  Power  spectral  density  of  Environment  I  at  SINR  =  -3  dB 
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In  addition,  the  TDOA  of  the  SOI  and  the  SNOI  are  separated  by  only  8  san^les, 
creating  some  spatial  overlap,  while  their  FDOA  vary  by  ^^256  separation  of 

the  parameters  of  these  two  signals  difficult  for  traditional  techniques  as  shown  in  Chapter 
n .  Note  as  a  further  example,  in  Figure  6-2  that  the  CAP  resolves  both  peaks  in  this  case. 
However,  the  SNOI  peak  is  dominant  resulting  in  an  incorrect  estimate  for  both  TDOA 
and  FDOA  for  the  SOI.  As  the  SINR  decreases,  the  SNOI  peak  dominates  more,  causing 
more  consistent  errors  in  the  estimate.  To  assess  the  performance  of  the  three  new 
cyclostationary  algorithms,  the  SNR  is  set  throughout  this  environment  at  0-dB  while  the 
SIR  starts  at  0-dB  and  decreases  through  -3  dB,  -10  dB  and  -13  dB,  representing  an 
increasingly  stronger  mterferer.  This  results  in  SINRs  that  cover  -  3  dB,  -  4.8  dB,  -10.4 
dB  and  - 13.2  dB. 

There  are  seven  collect  lengths  for  which  the  algorithms  are  tested.  Starting  with 
256  samples  which  represents  256  |i,s,  the  collect  lengths  increase  by  powers  of  2  to  16384 
samples,  representing  16.384  ms. 


Figure  6-2  -  CAF  of  Environment  I,  SINR  =  -3  dB  with  512  samples 
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3. 


Narrowband  Interference  Environment 


In  Environment  H,  the  interferer,  a  narrowband  SNOI,  has  10%  of  the  bandwidth 
of  the  SOL  In  this  case,  the  ability  of  the  algorithms  to  tolerate  more  focused  interference 
can  be  seen.  Rather  than  the  complete  spectral  and  temporal  overlap  of  the  vsddeband 
case,  which  results  in  some  of  the  SNOI  energy  lying  outside  the  bandwidth  of  the  SOI, 
the  energy  of  the  narrowband  SNOI  lies  completely  within  the  bandwidth  of  the  SOI.  This 
creates  a  different  type  of  problem  for  estimators  seeking  to  measure  the  parameters  of  a 
weak  SOI  in  the  presence  of  a  strong  much  more  spectrally  concentrated  SNOI.  In  the 
wideband  case,  some  of  the  energy  of  the  SNOI  can  be  eliminated  by  filtering  the  inputs  to 
the  algorithm  with  a  bandpass  filter  whose  center  frequency  is  equal  to  the  carrier  of  the 
SOI  and  with  bandwidth  equal  to  the  bandwidth  of  the  SOI.  In  the  narrowband  case,  this 
is  not  possible.  The  effects  of  that  energy  concentration  can  be  seen  in  Figures  6-3  and  6- 
4,  the  spectral  density  of  the  received  signal  and  an  example  plot  of  the  CAF  for  this 
environment. 

In  the  spectral  density  plot,  the  narrowband  SNOI  is  obvious  as  its  main  lobe  is 
much  narrower  and  rises  significantly  higher  than  the  overlapping  main  lobe  of  the  SOI. 
This  however  does  not  translate  well  to  a  separation  in  the  peaks  of  the  CAF.  In  fact,  the 
peaks  appear  to  combine  even  more  in  this  case  than  the  wideband  case.  At  this  point,  the 
SNOI  peak  is  again  dominant,  causing  an  incorrect  estimate  for  the  parameters  of  the  SOI. 
This  smearing  worsens  with  lower  SINR  and  eventually  causes  the  estimator  to  report  the 
presence  of  one  signal  instead  of  two.  Additionally,  the  estimates  often  obtained  are  good 
for  neither  the  SOI  nor  the  SNOI  due  to  the  combination  of  the  energy  from  each  of  the 
signals’  CAF  peaks. 

The  same  SNR,  SIR  and  SINR  levels  are  used  in  the  narrowband  case  as  were 
found  in  the  wideband  case.  This  allows  for  a  proper  comparison  of  the  performance  of 
the  algorithms  in  both  cases.  In  addition,  for  the  same  reason,  the  collect  lengths  remain 
the  same. 
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C.  TESTING 


1.  Stability  of  Estimates 

A  key  consideration  in  the  design  of  the  MCT  was  the  number  of  trials  per  case 
needed.  The  driving  factor  in  selecting  that  number  is  the  degree  of  stability  exhibited  by 
the  measure  of  the  variance  in  the  estimate  for  TDOA  and  FDOA,  produced  by  averaging 
the  error  results  over  all  the  trials  (see  next  section  for  discussion  of  averaging  technique). 
Initially,  400  trials  were  used  based  upon  testing  conducted  in  [22].  Unfortunately  and 
most  likely  due  to  the  combination  of  more  stressing  environments  and  different 
algorithms  tested,  the  400  trials  did  not  produce  a  stable  estimate  of  the  variance.  This  is 
obvious  in  the  inadequately  damped  sinusoid  of  Figure  6-5,  a  plot  of  the  cumulative  sum 
of  the  square  error  over  the  trial  number,  versus  the  trial  number.  An  adequate  number  of 
trials  would  produce  little  variation  in  the  final  few  values  on  the  plot  (Le.  a  well-damped 
sinusoidal  function).  For  these  algorithms,  400  trials  was  inadequate  for  such  a  curve. 


Figure  6-5  -  Stability  plot  -  TDOA,  400  trials,  Env.  I,  SINK  =  -3  dB,  512  samples 
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upon  further  investigation  of  [4],  it  appeared  that  1000  trials  produced  stable 
estimates  for  algorithms  similar  to  those  developed  here.  Therefore  with  an  additional  200 
trials  added  to  the  1000  for  margin,  the  resulting  stability  was  sufficient.  At  its  worst,  the 
variation  in  the  average  error  is  less  than  2%  for  the  TDOA  estimate  for  the  wideband 
interference  case,  SINK  =  -3  dB  and  512  samples.  At  its  best  it  is  less  than  0.25%  for  the 
estimate  of  FDOA  for  the  wideband  case,  SINR  =  -3  dB  and  1024  samples  and  several 
other  cases.  Both  stability  plots  for  these  extremes  can  be  seen  in  Figures  6-6  and  6-7 
respectively.  In  all  cases,  the  stability  of  the  FDOA  estimate  exceeded  that  of  the  TDOA 
estimate  by  a  factor  of  approximately  two. 


Figure  6-6  -  Stability  plot  -  TDOA,  1200  trials,  Env.  I,  SINK  =  -3  dB,  512  samples 
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Figure  6-7  -  Stability  plot  -  FDOA,  1200  trials,  Env.  I,  SINK  =  -3  dB,  1024  sauries 


2.  Products  of  Monte  Carlo  Testing 

Ultimately,  the  MCT  produces  a  measure  of  performance  for  each  collect  length  of 
each  of  the  eight  cases,  consisting  of  an  environment  and  a  SINK.  To  do  so,  each  trial  of 
the  MCT  produces  three  main  products:  a  D-fi  -  plane,  D  and  f<i,  the  MLEs  for  TDOA 
and  FDOA,  ^fd»  ^  TDOA  and  FDOA  estimates.  Each  of  the  1200 

trials,  that  make  up  each  case  produces  a  D-fd  -  plane,  similar  to  the  CAF  plane,  that  is  the 
manifestation  of  the  maximization  of  the  cyclostationary  TDOA-FDOA  algorithm  over 
time-delay,  D  and  Doppler,  fd.  Figure  6-8  shows  an  example  of  this  plane.  This  plane 
peaks  at  a  point  that  is  indicative  of  the  maximum  likelihood  estimates  for  the  TDOA  and 
FDOA  for  the  SOI. 
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Figure  6-8  -  Example  of  D-fd  -  plane  for  Env.  I,  SINR  =  -3  dB,  1024  samples 

The  maximum  likelihood  estimates  for  TDOA  and  FDOA  come  from  this  plane  as 
the  value  for  D  and  fd  that  produces  the  peak  in  the  D-fd  -  plane.  In  Figure  6-8,  the  MLE 
for  TDOA  can  be  seen  as  the  value  along  the  x-axis  that  corresponds  to  the  peak  in  the 
plane,  a  TDOA  of  47  samples  or  47  p.s  in  this  case.  The  MUE  for  FDOA  is  similarly 
found  as  the  value  on  the  y-axis  that  corresponds  to  the  peak  in  the  plane  - 19.53  kHz  for 
this  example.  Figures  6-9  and  6-10  show  the  slices  of  the  D_fi  -  plane  along  lines  of  time 
delay  and  frequency  shift  that  clearly  illustrate  the  TDOA  and  FDOA  respectively. 

Finally,  each  trial  produces  a  measurement  of  the  error  in  the  estimated  TDOA  and 
FDOA  by  finding  the  absolute  value  of  the  difference  in  the  measured  TDOA  and  FDOA 
versus  the  actual  TDOA  and  FDOA.  For  the  example  begun  in  Figure  6-8  and  continuing 
through  Figures  6-9  and  6-10,  those  errors  are  ^xdoa  =  3  samples  or  3  |j.s  and  fpDOA  =  0 
Kfertz.  Each  trial  thus  produces  a  D-fd  -  plane,  from  that,  an  MLE  for  both  TDOA  and 
FDOA  and  finally,  from  the  MI^s,  a  quantitative  measurement  of  the  error  m  the  estimate 
either  in  samples  for  TDOA  or  in  Hertz  for  FDOA. 
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Therefore,  for  each  Environment/SINR  case  at  each  of  the  seven  collect  lengths, 
1200  errors  for  TDOA  and  FDOA  estimation  are  produced  which  in  turn,  form  the  basis 
for  the  Root  Mean  Square  Error  (RMSE).  The  RMSE  is  the  ultimate  product  of  the  MCT 
and  is  the  parameter  used  to  compare  the  performance  of  the  newly  derived  algorithms  to 
that  of  the  CAP  and  the  stationary  CRB.  The  comparison  to  previously  derived 
cyclostationary  TDOA  algorithms  is  addressed  in  the  next  chapter  through  a  normalization 
of  the  RMSE. 

The  RMSE  for  TDOA  has  the  form  of  the  standard  Root  Mean  Square  Error 

where  L  is  the  number  of  trials,  D  is  the  actual  TDOA  in  seconds,  and  D^is  the  MLE  for 
the  TDOA  in  the  trial  in  seconds. 

The  performance  parameter  for  FDOA,  RMSEhxja  is  similar  to  RMSEtooa 
differing  only  in  the  units.  RMSEfdoa  has  the  form 


EMSEpDOA  d,B)  —  20  logjQ 


(6-2) 


where  again  L  is  the  number  of  trials,  is  the  actual  FDOA  in  Hertz,  and  f^^is  the  MLE  of 
the  FDOA  for  the  trial  in  Hertz. 

Thus,  with  two  Environments  and  four  SINR  levels  per  Environment  creating  eight 
distinctive  cases,  each  with  seven  collect  lengths  tested,  the  results  of  the  MCT  consist  of 
two  plots  for  each  case,  one  for  RMSEtdoa  and  one  for  RMSEftoa.  Each  of  these  sixteen 
plots  (eight  for  RMSEtdoa  and  eight  for  RMSEh^a)  contains  five  series  of  seven  points 
each.  The  five  series  represent  the  results  for  I2DMLE,  A2DMLE#1,  A2DMLE#2,  CAF 
and  stationary  CRB  with  the  seven  points  per  series  corresponding  to  the  results  obtained 
for  the  seven  collect  lengths  for  that  case. 
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D.  RESULTS 


1.  Performance  in  Environment  I  -  Wideband  Interference 

Figures  6-11, 6-12, 6-13  and  6-14  show  the  results  for  RMSEtooa  Environment  I, 
SINR  =  -3  dB,  SINR  =  -4.8  dB,  SINK  =  -10.4  dB  and  SINK  =  -13.2  dB  respectively. 
Following  that.  Figures  6-15, 6-16, 6-17  and  6-18  illustrate  the  RMSEfdoa  for 
Environment  I,  SINR  =  -3  dB,  SINR  =  -4.8  dB,  SINR  =  -10.4  dB  and  SINR  =  -13.2  dB 
respectively.  Note  that  the  x-axis  for  all  these  plots  is  the  collect  length  and  is  plotted  as  a 
log2-scaled  axis  to  simplify  the  visualization  of  data  taken  at  collect  lengths  that  range 
from  2®  to  2^*  samples. 


Figure  6-1 1  -  RMSEtdoa  -  Environment  I,  SINR  =  -3  dB 


RMSE(«,dB)  RMSE(t,dB) 


Figure  6-12  -  RMSEtooa  -  Environment  I,  SINK  =  -4.8  dB 


Figure  6-13  -  RMSEtdoa  -  Environment  I,  SINR  =  -10.4  dB 
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Figure  6-14  -  RMSEtdoa  -  Environment  I,  SINK  =  -13.2  dB 


Figure  6-18  -  RMSEfdoa  -  Environment  I,  SINR  =  -13.2  dB 

Because  the  algorithms  were  interpolated  on  the  frequency  axis  by  a  factor  of  8  at 
each  collect  length,  their  performance  is  tied  directly  to  the  minimum  frequency  resolution 
that  each  collect  length  employed  with  the  minimum  Epdoa  =  l*Af  per  trial  Thus,  the 
minimiim  erooA  ranges  from  l*Af = 488  Hz  at  AT  =  256  to  l*Af  =  7.6  Hz  at  iV  =  16384. 
The  level  of  interpolation  is  important  as  it  affects  the  RMSEfdoa  reported  in  the  next 
section.  With  more  interpolation  in  some  cases  it  is  possible  to  produce  smaller  errors  in 
the  estimate  of  FDOA  and  thus  lower  the  RMSE.  A  factor  of  8  was  chosen  as  a 
compromise  between  accuracy  and  computational  burden. 

2.  Performance  in  Environment  n  -  Narrowband  Interference 

Figures  6-19, 6-20, 6-21  and  6-22  show  the  results  for  RMSEtdoa  Environment  H, 
SINR  =  -3  dB,  SINR  =  -4.8  dB,  SINR  =  -10.4  dB  and  SINR  =  -13.2  dB  respectively. 
Finally,  Figures  6-23,  6-24, 6-25  and  6-26  illustrate  the  RMSEfdoa  for  Environment  H, 
SINR  =  -3  dB,  SINR  =  -4.8  dB,  SINR  =  -10.4  dB  and  SINR  =  -13.2  dB  respectively. 
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Figure  6-19  -  RMSEtdoa  -  Environment  H,  SINR  =  -3  dB 
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Figure  6-21  -  RMSEtdoa  -  Environment  II,  SINR  =  -10.4  dB 


Figure  6-22  -  RMSEtooa  -  Environment  H,  SINR  =  -13.2  dB 
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Figure  6-23  -  RMSEfdoa  -  Environment  H,  SINK  =  -3  dB 
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Figure  6-24  -  RMSEfdoa  -  Environment  H,  SINR  =  -4.8  dB 
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Vn.  OBSERVATIONS  AND  CONCLUSIONS 


A.  COMMENTS  ON  MONTE  CARLO  RESULTS 


1.  Ideal  vs.  Approxtinatioii 

I2DMLE  represents  the  theoretical  maximum  likelihood  algorithm  while 
A2DMLE#1  and  A2DM1JE#2  are  a  close  and  more  distant  approximation  respectively. 
Consequently,  I2DMLE  should  consistently  outperform  the  two  approximations,  which 
turns  out  to  be  the  case  in  all  instances.  The  main  objective  of  the  first  approximation, 
A2DMLE#1  was  to  maintain  all  six  terms  from  I2DMLE,  substitute  sin^)ly  for  the  pure 
signal  components  and  assess  the  performance  hoping  to  mimic  closely  that  of  the  ideal 
algorithm.  A2DMLE#1  does  in  fact  follow  closely,  the  trends  exhibited  by  the  ideal 
algorithm  hovering  generally  between  3  and  6  dB  below  the  performance  of  I2DMLE  in 
both  TDOA  and  FDOA.  The  two  obvious  exceptions  to  this  are  first,  in  the  higher  SINR 
cases,  the  slope  of  A2DMLE#1  RMSE  is  flatter  in  the  higher  collect  lengths  for  both 
TDOA  and  FDOA.  This  indicates  its  performance  is  still  improving  as  N  increases  but  to 
a  lesser  degree  than  I2DMLE.  Second,  in  the  narrowband  case,  A2DMLE#1  spears  to 
loose  more  performance  from  the  ideal  than  it  does  in  the  wideband  case  and  suffers  to  a 
greater  degree,  degradation  as  the  strength  of  the  narrowband  interference  increases  with 
increasing  SINR.  In  general  however,  it  is  a  very  good,  practical  approximation  to  the 
ideal 

A2DMLE#2  has  three  main  objectives  in  its  development:  (1)  keep  the  strongest 
terms  from  A2DMLE#1  to  maintain  as  much  as  possible  the  performance  of  A2DMLE#1, 
(2)  maintain  both  the  conjugate  and  non-conjugate  cyclostationary  exploitation  for 
flexibihty  and  (3)  minimize  the  computational  burden  by  limiting  the  number  of  terms  as 
much  as  possible.  With  these  dichotomous  goals  in  mind  the  resulting  three  terms  perform 
well  with  respect  to  A2DMLE#1  in  most  circumstances  and  reduce  the  computational 
burden  by  a  factor  of  2.  Its  performance  in  TDOA  is  more  even  than  either  the  ideal  or 


107 


the  first  approximation  progressing  with  a  near  linear  slope  on  the  RMSE  curves.  It  too 
slows  its  increase  in  performance  as  the  collect  length  increases  to  16384  samples  as  is 
evident  by  the  increase  in  the  slope  of  the  performance  line  at  the  longest  coDect  lengths. 
In  the  narrowband  case  for  TDOA,  it  is  surprisingly  well  behaved,  with  performance  near 
or  in  excess  of  A2DMLE#1.  Finally,  in  FDOA,  the  same  trends  found  in  TDOA  continue 
but  with  an  apparent  threshold  around  SINR  =  -13  dB  where  it  again  rises  back  above  the 
A2DMLE#1  line.  In  general,  A2DMLE#2  is  another  good  approximation  resulting  in 
reduced  computational  burden  over  A2DMLE#1  and  appears  to  endure  narrowband 
interference  more  successfully. 

2.  New  Cyclostationary  Estimators  and  the  Cross  Ambiguity  Function 

Comparison  with  the  CAF  illustrates  the  value  of  cyclostationary  techniques  in 
these  heavy  interference  environments.  In  general  for  TDOA,  at  shorter  collect  lengths 
and  higher  SINR,  the  CAF  outperforms  all  but  the  ideal  algorithm  and  in  some  cases, 
outperforms  all  three.  But  as  collect  length  increases,  all  three  cyclostationary  algorithms 
begin  to  outperform  the  CAF  until  at  the  longest  collect  lengths,  the  difference  in 
performance  is  as  great  as  15  to  16  dB  (seconds)  between  I2DMIJE  and  the  CAF  at  SINR 
=  -13  dB.  Also,  as  SINR  decreases,  the  CAF  rapidly  falls  off,  until  at  SINR  =  -13  dB,  all 
three  cyclostationary  algorithms  outperform  the  CAF  at  every  collect  length.  There 
seemed  a  slight  difference  in  the  CAF  performance  between  the  wideband  and  the 
narrowband  environments  with  the  performance  worse  with  narrowband  interference. 

For  FDOA,  the  CAF  actually  performs  better  than  all  three  cyclostationary 
algorithms  for  a  longer  period  than  was  true  for  TDOA.  At  higher  SINR  in  the  wideband 
case,  the  CAF  outperforms  all  three  cyclostationary  algorithms  until  I2DMLE  crosses  at  N 
=  2048.  Following  that,  the  two  approximations  cross  at  AT  =  4096  and  all  three  continue 
to  outperform  the  CAF  from  that  point.  At  higher  SINR  in  the  narrowband  case,  the  CAF 
maintains  superiority  through  N  =  4096  when  all  three  cyclostationary  algorithms  exceed  it 
performance.  Finally,  for  lower  SINR  in  both  the  wideband  and  narrowband 
environments,  the  CAF  performance  deteriorates  quickly  to  the  point  at  which  its  results 
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are  consistently  in  error  and  dominated  completely  by  the  parameters  of  the  interfering 
signal  This  is  dear  in  the  SINK  =  -13  dB  cases  where  the  CAP  is  significantly  worse  than 
all  three  cyclostationary  algorithms. 

3.  New  Cydostationary  Estimators  and  the  Cramer-Rao  Bounds 

The  three  new  cydostationary  algorithms  perform  well  in  con^arison  to  the  lower 
bound  predicted  by  the  stationary  CRB.  Because  it  is  derived  under  the  assumption  that 
the  environment  is  pure  Gaussian  and  stationary  in  nature  and  because  cydostationary 
theory  actually  requires  non-stationary  signals  in  order  to  achieve  the  signal-selectivity  and 
ultimately  the  processing  gains  needed,  the  stationary  CRB  is  only  useful  as  a  qualitative 
conparison.  While  it  will  not  produce  empirical  evidence  as  to  the  theoretical  maximum 
contribution  that  cydostationary  techniques  make  over  stationary,  this  comparison  does 
produce  a  notion  of  the  power  of  these  new  algorithms  and  the  contribution  they  can  make 
in  practice. 

In  TDOA,  the  power  of  the  cydostationary  techniques  is  not  as  evident  in  the 
shorter  collect  lengths  and  higher  SINR  levels  as  it  is  in  the  longer  collect  lengths  and  low 
SINR  levels.  The  RMSE  of  the  three  algorithms  all  exceed  that  of  the  CRB  for  SINR  =  -3 
dB,  SINR  =  -4.8  dB  and  SINR  =  -10.4  dB  for  collect  lengths  under  4096  samples. 
I2DMLE  begins  to  perform  better  than  CRB  predictions  at  8192  samples  for  these  three 
SINR  levels.  The  same  is  true  at  the  shorter  collect  length  of 4096  samples  for  SINR  =  - 
10.4  and  SINR  =  -13.2.  Neither  approximation  exceeds  the  performance  of  the  CRB 
predictions  however,  A2DMLE#1  does  approach  the  CRB  as  the  collect  length  increases 
to  16384  samples.  Both  wideband  interference  and  narrowband  interference  produce 
similar  results. 

In  FDOA,  the  reverse  pattern  appears  in  the  wideband  interference  environment. 
I2DMLE  performs  exceptionally  well  at  SINR  =  -3  dB  exceeding  the  bound  at  4096 
samples  and  continuing  out  to  16384  samples.  However,  as  the  SINR  drops,  so  does  the 
performance  of  I2DMLE  with  respect  to  the  CRB.  By  SINR  =  -13  dB,'  only  at  W  =  16384 
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does  I2DMLE  exceed  the  bound.  The  two  approximations  again  never  exceed  the  bound 
however,  both  appear  to  approach  the  CRB  as  N  increases  to  16384  and  beyond. 

In  the  narrowband  case,  the  relative  performance  of  I2DMLE  is  worse  yet.  At 
SINR  =  -3  dB,  only  N  =  16384  produces  a  result  under  the  bound.  This  is  true  for  each  of 
the  SINR  levels  with  the  difference  between  I2DMLE  and  the  CRB  narrowing  as  SINR 
increases.  This  leads  naturally  to  the  conclusion  that  the  cyclostationary  algorithms  suffer 
more  in  narrowband  interference  than  in  wideband  interference,  a  trait  they  share  in 
conunon  with  the  CAP  according  to  Stein  in  [7]. 

The  effects  of  interpolation  on  the  frequency  axis  definitely  affect  the  performance 
of  the  algoiithm  Because  only  a  factor  of  8  was  used  to  interpolate  at  each  collect  length, 
the  minimum  error  possible  for  each  MLEfdoa  range  from  the  maximum  of 488  Hz  at  iV  = 
256  to  the  minimum  of  7.6  Hz  at  =  16384.  It  is  likely  that  the  algorithms,  the  CAP 
included,  would  produce  RMSE  much  lower  given  a  greater  interpolation  factor  and  thus 
approach  or  exceed  the  CRB  in  more  cases  (with  the  exception  of  the  CAP  which  of 
course  cannot  exceed  the  CRB).  This  increased  interpolation,  which  increases  the 
con^utational  burden  as  the  square  of  the  increase  in  the  interpolation  factor,  was  not 
pursued.  It  should  simply  be  noted  that  the  RMSE  for  three  cyclostationary  algorithms 
and  the  CAP  are  probably  conservative  estimates  in  light  of  this  conclusion. 

4.  Normalization  and  Comparison  to  Herit^e  Algorithms 

In  order  to  compare  the  RMSE  results  from  above  to  those  for  previously  derived 
cyclostationary  TDOA  algorithms  contained  in  [4]  and  [22],  it  is  necessary  to  normalize 
the  RMSE  in  some  way,  in  order  to  account  for  the  difference  in  the  bandwidth  of  the  SOI 
tested  here  versus  that  used  in  each  of  the  two  heritage  developments.  The  simplest 
method  for  doing  involves  adjusting  the  RMSE  results  for  the  difference  in  bandvridth. 

The  bandwidth  of  the  SOI  used  to  test  the  three  new  cyclostationary  algorithms  above  is 
187.5  kHz.  To  calculate  the  performance  in  the  same  environments  for  a  SOI  with 
bandwidth  of  1.875  kHz,  a  simple  adjustment  to  the  RMSE  results  is  all  that  is  needed. 
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Equation  (7-1)  shows  that  the  single  correction  factor  to  accommodate  SOI  of  different 
bandwidths. 
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The  net  effect  is  a  decrease  in  the  RMSE  for  SOI  of  greater  bandwidth  and  an  increase  for 
SOI  of  lesser  bandwidth.  For  example,  for  a  SOI  of  1.875  MHz,  the  RMSE  would 
decrease  by  -20  dB.  For  a  SOI  with  bandwidth  18.75  kHz,  the  increase  is  20  dB. 

With  this  is  mmd,  adjustment  can  be  made  for  the  difference  in  the  bandwidth  of 
the  SOI  here  and  that  contained  in  [22].  Although  the  environments  and  SINR  levels 
remain  different,  a  qualitative  assessnoent  of  the  comparative  performance  is  possible. 

Adjusting  the  results  of  Chapter  VI  for  the  difference  in  bandwidth,  an  adjustment 
of +3.5  dB,  it  is  clear  that  while  comparable,  the  performance  of  I2DMLE  and  the  two 
approximations,  A2DMLE#1  and  A2DMLE#2,  lag  behind  the  similar  heritage  algorithms 
by  2  to  4  dB  in  similar  SINK  conditions.  At  &st  glance  it  would  appear  that  the  heritage 
algorithms  outperform  those  newly  derived,  but  at  closer  examination,  the  interference 
environments  and  specifically  the  type  of  interference  used  here  clearly  are  more  difficult 
an  environment  in  which  to  operate.  Thus,  the  new  algorithms  more  than  likely  perform 
nearly  equal  to  those  derived  and  tested  in  [22], 

Note  that  these  comparisons  hold  true  only  for  TDOA.  While  some  work  into  the 
tolerance  of  these  the  TDOA  algorithms  in  [22]  to  various  FDOA  has  been  reported,  little 
development  and  even  less  testing  has  occurred  in  true  cyclostationary  TDOA-FDOA, 
leaving  no  viable  conparison  to  be  made  for  the  FDOA  results. 


5.  Anollier  Look  and  Perfonnance  vs.  Collect  Length 

Another  method  for  evaluating  the  performance  as  collect  length  increases  comes 
firom  examining  the  slope  of  the  lines  plotted  as  RMSE  *  collect  length  vs.  collect  length. 
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From  this,  it  is  easy  to  tell  if  the  performance  of  an  algorithm  is  increasing  more  or  less 
than  the  collect  length  is  increasing  which  in  turn,  gives  some  sense  of  the  added  benefit  of 
collecting  and  processing  more  signal  data  at  the  cost  of  computational  burden.  If  an 
algorithm  has  a  TDOA  performance  line  with  shallow  slope,  more  data  improves 
performance  at  less  cost  than  an  algorithm  with  a  relatively  steep  slope.  For  FDOA,  as 
explained  below,  algorithms  with  slopes  equal  to  or  less  than  the  slope  of  the  CRB  line  are 
increasing  performance  better  than  or  proportional  to  the  increase  in  collect  length.  Those 
with  greater  slopes  are  not. 

An  abbreviated  examination  of  this  phenomena  appears  in  Figures  7-1  through  7-8, 
the  RMSE  ♦  collect  length  vs.  collect  length  plots  for  TDOA  and  FDOA  estimation  in 
Environments  I  and  n  at  SINK  =  -  3  dB  and  -13  dB.  These  eight  plots  represent  the 
extremes  of  the  performance  investigated  here  and  thus,  conclusions  about  the  trends  in 
the  performance  of  the  algorithms  can  be  noted  without  the  need  for  plotting  the 
remaining  eight  middle  SINR  cases.  Again  the  x-axis  containing  the  collect  length  is 
plotted  as  a  log2-scaled  axis. 

Note  also  that  because  the  relationships  between  collect  length,  the  minimum 
£fdoa»  find  RMSEfdoa  fire  more  complex  for  FDOA  than  for  TDOA,  as  explained  above, 
the  RMSEhx)a  plots  below  are  plotted  as  logio  vs.  log2  plots  to  obtain  a  more  linear 
relationship.  Again,  while  a  slope  of  zero  does  not  indicate  performance  increases  equally 
with  collect  length  for  these  FDOA  plots,  it  is  obvious  from  the  slope  of  the  CRB  lines 
what  constitutes  a  consistent  increase  in  performance  with  increasing  collect  length.  The 
CRB  lines  for  FDOA  are  thus  the  benchmark  here. 

The  results  in  Figures  7-1  through  7-8  seem  to  indicate  that  the  cyclostationary 
algorithms  have  minor  thresholds  at  the  higher  SINR  level  This  is  displayed  in  the  stair¬ 
step  pattern  seen  in  both  Figures  7-1  and  7-5.  This  pattern  is  also  evident  to  a  tesser 
extent  in  the  RMSEhx)a  plots  in  Figures  7-2  and  7-6  also  confirming  that  the  performance 
in  time  and  frequency  measurement  is  closely  related. 

For  SINR  =  -3  dB,  with  the  exception  of  A2DMLE#2  which  must  incur  a 
significant  increase  in  computational  burden  in  order  to  significantly  increase  its 
performance,  indicated  by  the  sleep  slope  of  its  line  after  N  =  8192,  the  cyclostationary 
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algorithms  appear  to  mimic  the  CRB.  Both  I2DMLE  and  A2DMLE#1  foflow  the  slope  of 
the  CRB,  concluding  at  iV  =  16384  with  roughly  the  same  slope  as  the  CRB.  Two 
conclusions  may  be  drawn  from  this:  (1)  I2DMLE  and  A2DMLE#1  will  gain  little  in 
performance  relative  to  the  CRB  with  increasing  collect  length  and  (2)  these  two 
cyclostationary  algorithms  require  the  same  increase  in  computational  burden  in  order  to 
increase  performance  by  a  set  amount. 

This  stair-step  patterns  of  thresholds  in  both  TDOA  and  FDOA  are  smoothed  in 
the  low  SINR  level  plots  of  Figures  7-3, 7-4, 7-7  and  7-8,  dampened  by  the  significant 
increase  in  interference  at  this  level  For  this  SINR  level  in  TDOA  performance,  the  three 
cyclostationary  algorithms  again  follow  the  slope  of  the  CRB  with  the  exception  of  the 
two  approximations  a.tN=  16384  in  the  narrowband  interference  case.  In  FDOA,  the 
three  cyclostationary  algorithms  actually  conclude  with  steeper  negative  slopes  than  the 
CRB.  Thus  though  they  rarely  exceed  the  bound  up  to  N  =16384,  it  appears  that  for 
longer  collect  lengths,  all  three  would  eventually  exceed  the  CRB.  In  addition,  it  would 
appear  that  increasing  collect  length  improves  the  ability  of  the  cyclostationary  algorithms 
to  measure  FDOA  more  so  than  the  improvement  that  could  be  expected  from  a 
stationary-based  estimator  bounded  by  the  CRB. 


Figure  7-1  -  Performance  as  N  increases  -  TDOA,  Environment  I,  SINR  =  -3 
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Figure  7-4  -  Performance  as  N  increases  -  FDOA,  Environment  I,  SINK  =  -13 
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Figure  7-5  -  Performance  as  N  increases  -  TDOA,  Environment  n,  SINR  =  -3 
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Figure  7-8  -  Performance  as  N  increases  -  FDOA,  Environment  H,  SINK  =  -13 


6.  Perfonnance  as  a  Function  of  Signal  to  Interference  and  Noise  Ratio 

Another  comparison  that  is  particularly  constructive  is  that  of  performance  as  a 
function  of  SINK.  One  of  the  theoretical  advantages  of  cyclostationary  techniques  is  the 
highly  signal-selective  properties  offered  by  algorithms  that  take  advantage  SOI-unique 
cycle  frequencies  in  order  to  eliminate  energy  in  the  processing  bandwidth  that  is  not 
attributable  to  the  SOL  Thus,  cyclostationary  algorithms  should  degrade  in  performance 
as  SINK  decreases  much  more  gracefully  than  traditional  techniques  and  the  stationary 
CRB. 

The  fourteen  plots  in  Figures  7-9  through  7-22  illustrate  the  RMSEtooa  as  a 
function  of  SINR  for  each  of  the  seven  collect  lengths.  Figures  7-9  through  7-15  are  for 
Environment  I  and  Figures  7-16  thought  7-22  correspond  to  Environment  n. 
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Figure  7-9  -  RMSEtdoa  vs.  SINR,  Environment  I,  =  256 


Figure  7-10  -  RMSEtdoa  vs.  SINR,  Environment  I,  N  =  512 
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Figure  7-1 1  -  RMSEtdoa  vs.  SINR,  Environment  I,  N=  1024 


Figure  7-12  -  RMSEtdoa  vs.  SINR,  Environment  I,  N  =  2048 


Figure  7-13  -  RMSEtdoa  vs.  SINR,  Environment  I,  =  4096 


Figure  7-14  -  RMSEtooa  vs.  SINR,  Environment  1,N=  8192 
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Figure  7-15  -  RMSEtooa  vs.  SINK,  Environment  1,N=  16384 


Figure  7-16  -  RMSEtdoa  vs.  SINR,  Environment  n,  =  256 
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Figure  7-17  -  RMSEtdoa  vs.  SINK,  Environment  H,  =  512 


Figure  7-18  -  RMSEtdoa  vs.  SINR,  Environment  IL,N=  1024 
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Figure  7-19  -  RMSEtdoa  vs.  SINK,  Environment  n,  iV = 2048 


Figure  7-20  -  RMSEidoa  vs.  SINR,  Environment  II,  iV  =  4096 
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Figure  7-21  -  RMSEtdoa  vs.  SINK,  Environment  H,  =  8192 


Figure  7-22  -  RMSEtdoa  vs.  SINR,  Environment  H,  IV  =  16384 
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The  most  obvious  result  illustrated  by  these  plots  is  the  confirmation  of  the 
conclusion  from  above  with  respect  to  the  cyclostationary  algorithms’  performance  in 
shorter  collect  lengths  and  lower  SINK  levels.  Beginning  with  Figures  7-9  and  7-16,  for  N 
=  256,  the  cyclostationary  algorithms  perform  well  above  the  CRB,  in  the  same  range  as 
the  CAF  and  do  not  increase  their  performance  relative  to  these  two  stationary 
benchmarks  as  SINR  decreases.  However,  as  is  clear  in  Figures  7-13  and  7-20,  by  = 
4096,  the  true  contribution  of  the  cyclostationary  algorithms  becomes  evident.  As  SINR 
deteriorates,  these  three  algorithms  continue  to  perform  well  and  significantly  increase 
their  performance  relative  to  the  stationary  techniques.  By  =  16384,  in  Figures  7-15 
and  7-22,  the  cyclostationary  algorithms  significantly  outperform  the  CAF  and  appear  to 
suffer  less  degradation  in  performance  as  SINR  decreases  than  even  the  CRB. 

Following  the  same  theme  for  FDOA,  it  too  is  instructive  to  understand  the 
sensitivity  of  the  algorithms  to  the  SINR  level  at  each  of  the  collect  lengths.  Figures  7-23 
through  7-29  show  the  RMSEfdoa  vs.  SINR  for  Environment  I  and  Figures  7-30  through 
7-36  illustrate  the  results  for  Environment  II. 


Figure  7-23  -  RMSEfdoa  vs.  SINR,  Environment  1,N=  256 


125 


nil8E(«.dB) 


100.0 


70.0-1 


80.0 


•3.0 


— 5K — Cramer-Rao 

•10.4  -13^2 


8INR  (dB) 


Figure  7-24  -  RMSEh)oa  vs.  SINR,  Environment  I,  iV  =  512 


Figure  7-25  -  RMSEfdoa  vs.  SINR,  Environment  I,  =  1024 


126 


nusEd^dB) 


Figure  7-26  -  RMSEedoa  vs.  SINR,  Environment  I,  iV=  2048 


Figure  7-27  -  RMSEfdoa  vs.  SINR,  Environment  I,  N  =  4096 
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Figure  7-28  -  RMSEfdoa  vs.  SINR,  Environment  1,N=  8192 


Figure  7-29  -  RMSEfdoa  vs.  SINR,  Environment  I,  iV  =  16384 
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Figure  7-30  -  RMSEfdoa  vs.  SINK,  Environment  II,  =  256 


Figure  7-31  -  RMSEfdoa  vs.  SINR,  Environment  H,  iV  =  512 
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Figure  7-32  -  RMSEfdoa  vs.  SINR,  Environment  H,  iV  =  1024 


Figure  7-33  -  RMSEhxja  vs.  SINR,  Environment  11,N  =  2048 
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Figure  7-34  -  RMSEfdoa  vs.  SINR,  Environment  I1,N  =  4096 


Figure  7-35  -  RMSEfdoa  vs.  SINR,  Environment  TL,N=  8192 
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Figure  7-36  -  RMSEhxja  vs.  SINK,  Environment  11,N=  16384 

From  these  RMSEfdoa  vs.  SINK  plots,  several  conclusions  may  be  drawn.  First, 
the  gap  between  the  implemented  algorithms  and  the  CRB  at  short  collect  lengths  is  very 
obvious  in  Figures  7-23  and  7-30,  the  gap  a  bit  larger  in  7-30  due  to  the  more  difficult 
narrowband  interference  environment  as  noted  above.  In  addition,  as  SINR  decreases  at 
this  short  collect  length,  the  implemented  algorithms  degrade  more  rapidly  than  the  CRB. 
Through  the  next  sets  of  plots.  Figures  7-24  through  7-28  and  7-31  through  7-35,  the 
gradual  progression  in  the  performance  of  the  cyclostationary  algorithms  can  be  seen.  By 
N  =  8192  in  Figures  7-28  and  7-35,  the  cyclostationary  algorithms  have  begun  to  approach 
the  CRB,  now  degrading  at  the  same  rate  or  even  better  than  the  CRB.  In  this  same  plot, 
the  fate  of  the  CAF  is  obvious  with  its  RMSE  well  above  those  of  the  cyclostationary 
algorithms. 

Finally,  in  Figures  7-29  and  7-36,  the  cyclostationary  algorithms  all  conclude  at 
SINR  =  -13  with  less  degradation  than  the  CRB  indicating  that  as  SINR  decreases  further, 
all  three  algorithms  will  soon  exceed  the  lower  bound  set  by  the  CRB. 
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7.  Siuninaiy 


Conclusions  in  six  major  areas  are  available  as  a  result  of  the  MCT  presented  in 
Chapter  VI.  These  areas  include  the  effects  of  approximation  on  the  ideal  cyclostationary 
two-dimensional  maximum  likelihood  algorithm,  I2DMLE;  a  comparison  of  the  new 
cyclostationary  algorithms  to  their  stationary  counterpart,  the  CAF;  a  comparison  of  the 
new  cyclostationary  algorithms  to  the  stationary  theoretical  lower  bound  on  the  accuracy 
of  TDOA  and  FDOA  estimation,  the  CRB;  a  normalization  and  con^arison  of  the  newly 
derived  algorithms  to  heritage  cyclostationary  TDOA  algorithms;  a  different  illustration  of 
performance  as  collect  length  increases;  and  an  investigation  of  performance  as  SINR 
decreases. 

The  effects  of  approximation  on  I2DMLE  can  be  seen  in  the  performance  of 
A2DMLE#1  and  A2DMI^#2.  The  first  approximation,  A2DMLE#1,  followed  closely 
the  performance  of  the  ideal  algorithm,  differing  in  RMSEtdoa  and  RMSEfdoa  by  2  to  4 
dB  across  nearly  all  cases.  It  is  deemed  a  good,  practical  approximate  cyclostationary  ML 
TDOA-FDOA  algorithm  due  to  its  reasonable  performance  relative  to  the  ideal  and  its 
lack  of  dependence  on  a  priori  information  like  the  spectral  correlation  fiinction  of  the 
SOI.  A2DMLE#2  is  a  compromise  between  performance  and  computational  conq)lexity. 
Its  performance  is  3  to  5  dB  above  that  of  A2DMLE#1  and  it  lacks  the  steep  sloping 
performance  vs.  collect  length  at  the  longer  collect  lengths.  However,  in  its 
implementation,  it  enjoys  computational  savings  of  nearly  2.  Thus,  it  is  a  good 
compromise  between  performance  and  computational  burden. 

An  three  cyclostationary  algorithms  performed  well  in  comparison  to  the  CAF.  In 
TDOA,  the  CAF  outperforms  the  cyclostationary  implementations  at  very  short  collect 
lengths  of  256  or  512  samples,  but  only  at  higher  SINR.  When  the  signal  environment 
deteriorates,  the  cyclostationary  algorithms  continue  to  improve  their  performance  as 
collect  length  increased  while  the  performance  curve  of  the  CAF  flattens  out  as  a  result  of 
the  increased  interference.  Thus,  at  low  SINR  and  longer  collect  lengths,  the  three 
cyclostationary  algorithms  outperformed  the  CAF  by  as  much  as  16  dB-(seconds).  In 
FDOA,  the  cyclostationary  algorithms  are  outperformed  at  higher  SINR  by  the  CAF  until 


N  =  4096,  quite  a  bit  longer  than  in  TDOA.  However  again  as  the  signal  environment 
deteriorates,  the  performance  of  the  cyclostationary  algorithms  does  not  deteriorate  as 
rapidly  as  that  of  the  CAF  resulting  in  a  42  dB  (Hertz)  gap  at  SINR  =  -13  dB  and  = 
16384.  Clearly,  the  advantage  of  the  cyclostationary  algorithms  over  the  CAF  in  heavy 
interference  environments  is  evident  by  these  results. 

Performance  of  the  cyclostationary  algorithms  in  comparison  to  the  stationary 
theoretical  lower  bounds  on  TDOA  and  FDOA  estimation  set  by  the  CRB  illustrate  the 
contribution  that  the  signal-selectivity  and  interference  tolerance  features  of  the 
cyclostationary-based  algorithms  make.  In  TDOA,  I2DMLE  exceeds  the  CRB  for  larger 
N  and  the  two  approximations  approach  the  CRB  as  N  approaches  16384  samples.  This 
indicates  that  with  increasing  N,  all  three  cyclostationary  algorithms  exceed  the  CRB.  In 
FDOA,  again,  I2DMLE  exceeds  the  CRB  for  N  >  4096  and  SINR  =  -3  dB.  As  SINR 
decreases,  however,  the  number  of  samples  required  for  the  cyclostationary  algorithms  to 
exceed  the  CRB  increases.  At  SINR  =  -13  dB,  only  I2DMLE  exceeds  the  CRB  and  only 
with  N  =  16384.  But,  as  in  TDOA,  all  three  approach  the  CRB  as  N  approaches  16384 
samples.  This  clearly  indicates  that  in  FDOA  as  in  TDOA,  as  N  exceeds  16384,  each  of 
the  cyclostationary  algorithms  exceeds  the  CRB. 

Normalization  of  these  results  for  SOI  of  different  bandwidths  is  accomplished 
with  a  simple  adjustment.  This  adjustment  amounts  to  -i-3.5  dB  for  comparison  to 
algorithms  outlined  in  [22].  With  this  adjustment,  the  three  new  cyclostationary 
algorithms  performed  within  3  to  5  dB  of  the  heritage  algorithms.  The  more  challenging 
interference  environments  used  m  the  MCT  above  can  account  for  at  least  some  of  that 
performance  gap.  Thus,  the  results  obtained  here  are  consistent  with  those  achieved  with 
previously  derived  estimators. 

Performance  as  a  function  of  collect  length  is  an  important  parameter,  particularly 
in  making  implementation  decisions  like  increasing  the  number  of  samples  processed.  By 
plotting  the  RMSE  *  collect  length  vs.  collect  length,  it  is  readily  apparent  which 
algorithms  benefit  most  from  increasing  collect  lengths  to  be  processed.  With  the 
exception  of  A2DMLE#2  in  TDOA  estimation,  the  cyclostationary  algorithms  equaled  or 
exceeded  any  performance  gain  made  by  the  CRB  due  to  an  increase  in  collect  length. 
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Thus,  for  improved  performance,  at  the  expense  of  computation  time,  all  three 
cyclostationary  algorithms  improve  nearly  linearly  with  increased  collect  length  with  the 
exception  noted  above. 

Finally,  looking  at  performance  as  a  function  of  SINR  again  illustrates  the  signal- 
selectivity  and  interference  tolerant  characteristics  of  the  cyclostationary  algorithms. 

While  the  CAF  and  it’s  theoretical  lower  bound  on  performance  set  by  the  CRB  degrade 
significantly  as  the  SINR  decreases  to  -13  dB,  the  cyclostationary  algorithms  do  not  suffer 
nearly  the  same  losses  in  performance.  In  fact,  all  three  cyclostationary  algorithms 
{^proach  or  exceed  the  CRB  as  SINR  decreases  to  —13  dB.  Again,  another  indication  of 
the  power  that  these  cyclostationary  algorithms  bring  to  highly  corrupt  signal 
environments. 

The  results  of  this  work  clearly  illustrate  three  joint  TDOA-FDOA  estimators, 
optimized  in  a  maximum-likelihood  sense,  that  take  full  advantage  of  cyclostationary 
signal-selectivity  and  interference  tolerance.  They  perform  well  in  con^)arison  to  both  the 
CAF  and  the  stationary  CRB,  indicating  the  potential  performance  gains  cyclostationary 
techniques  have  to  offer  over  traditional  techniques.  In  addition,  their  performance  in 
TDOA  is  comparable  to  heritage  TDOA  algorithms  when  adjusted  for  SOI  bandwidth. 
Finally,  they  continue  to  enjoy  performance  gains  as  N  increases  and  also  as  SINR 
decreases,  not  normally  seen  in  traditional  techniques. 

B.  ORIGINAL  CONTRIBUTIONS 

The  original  contribution  of  this  research  consists  of  the  new  three  algorithms 
developed  and  the  performance  testing  of  these  new  algorithms.  I2DMLE  represents  the 
first  development  of  a  joint  maximum-likelihood  TDOA-FDOA  estimator  for  digital 
cyclostationary  communications  signals.  It  is  derived  from  a  modified  version  of  Model  C, 
the  simple  Doppler  model,  which  has  yet  to  be  used  as  the  basis  for  the  development  of 
any  cyclostationary  TDOA  or  TDOA-FDOA  estimators.  As  noted  above,  it  is  an  ideal 
algorithm  requiring  the  a  priori  knowledge  of  the  spectral  correlation  function  of  the  SOI 
in  order  to  be  implemented.  This  situation  being  difficult  to  accommodate  in  practice,  a 


substitution  scheme  is  necessary  to  form  approximations.  Whfle  an  infinite  number  of 
substitutions  is  available  to  produce  an  infinite  number  of  approximations,  A2DMLE#1 
and  A2DMLE#2  represent  two  original  approximations  to  I2DMLE.  Thus,  they  are  the 
first  two  practical  cyclostationary  joint  ML  estimators  for  TDOA  and  FDOA  requiring 
only  a  cycle  frequency  from  the  SOI  rather  than  the  entire  spectral  correlation  function. 

The  con]|>arison  of  the  performance  of  these  three  original  algorithms  to  that  of  the 
CAP  and  the  stationary  CRB  predictions  also  marks  the  first  time  any  cyclostationary 
TDOA-FDOA  estimation  algorithm  has  been  evaluated  m  such  detail  against  traditional 
techniques.  While  extensive  work  has  been  done  in  the  area  of  cyclostationary  TDOA 
estimators,  the  work  in  cyclostationary  TDOA-FDOA  determination  has  been  restricted  to 
the  investigation  of  the  sensitivity  of  cyclostationary  TDOA  estimators  to  FDOA  and  to 
the  modification  of  TDOA  algorithms,  derived  fi*om  Model  A,  to  accommodate  FDOA. 
Thus,  the  performance  of  I2DMLE,  A2DMLE#1  and  A2DMLE#2  is  the  first  benchmark 
to  which  future  cyclostationary  TDOA-FDOA  algorithms  can  be  compared. 


C.  AREAS  FOR  FURTHER  INVESTIGATION 


1.  Hme-varying  Doppler 

The  effects  of  the  more-realistic  case  of  time- varying  Doppler  pose  both  a 
challenging  and  intriguing  puzzle.  Its  first  influence  in  the  TDOA/FDOA  problem  comes 
in  the  form  of  various  unknown  effects  on  the  calculation  of  the  SCF.  How  time- varying 
Doppler  changes  the  SCF  surface  and  affects  the  determination  of  cycle-frequencies-of- 
interest  has  been  investigated  little  if  any.  Beyond  that,  its  effect  on  the  two-dimensional 
optimization  problem  derived  and  tested  here  becomes  even  less  clear.  Certainly, 
allowances  have  to  be  made  and  a  proper  choice  for  determining  the  instantaneous  FDOA 
would  have  to  be  determined.  In  addition,  other  more  subtle  changes  certainly  will  be 
discovered  once  a  detailed  analysis  is  made. 
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2. 


Various  Spectral  Correlation  Function  Substitution  Options 


As  noted  in  Chapter  V,  the  derivation  concludes  with  the  logical  selection  of 
substitutions  for  the  SCF  of  the  SOI  based  upon  the  two  observation  variables.  Those 
substitutions,  being  driven  primarily  by  arbitrarily  chosen  rules  could  be  quite  different 
should  another  set  of  rules  be  selected  based  on  one  of  the  other  possible  motivations 
outlined  in  Chapter  V  as  well  These  substitutions  could  change  the  performance  of  the 
algorithm  significantly  in  addition  to  affecting  the  digital  implementation  aspects  of  the 
problem.  Both  accurate  determination  of  the  TDOA  and  FDOA  in  addition  to  tractable 
implementation  on  general  purpose  or  DSP  hardware  is  relevant  and  deserves  additional 
attention.  Rnally,  the  inclusion  of  the  phase  and  anq)litude  mismatch,  for  highly  disparate 
and  non-coherent  platforms  does  add  the  additional  problem  of  estimating  the  mismatches 
first  (see  below).  With  this  in  mind,  the  substitutions  could  be  made  to  simplify  this 
process  or  minimize  the  reliance  on  these  parameters  that  themselves  may  have  suboptimal 
estimates. 


3.  Amplitude  and  niase  Mismatch 

While  arguments  outlined  in  Chapter  V  supported  the  disregard  of  the  mismatch  in 
amplitude  and  phase  between  the  receivers,  there  may  exist  an  instance  where  the 
coherence  between  platforms  is  poor  enough  to  play  a  significant  role  in  the  performance 
of  the  algorithm.  While  some  work  on  specific  algorithms  has  been  done  (SSPI  contract), 
a  study  that  concentrates  specifically  on  the  general  effects  of  amplitude  and  phase 
mismatch  at  RF  fi*equencies  on  TDOA  and  FDOA  determination  using  cyclostationary 
techniques  has  not  been  published.  This  work  could  ideally  lay  the  rigorous  mathematical 
foundation  for  the  inclusion  or  exclusion  of  the  mismatch  in  realistic  signal  models,  thus 
making  a  fundamental  contribution  to  this  field. 
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4. 


Modulation  Types 


As  the  spectram  crowds  evermore,  many  new  and  more  complex  modulation  types 
are  in  use.  The  work  presented  here  tested  only  phase-shift  keying  schemes,  specifically, 
BPSK  with  the  logical  and  mathematical  extension  to  QPSK.  Higher-order  PSK  schemes 
as  well  as  the  myriad  of  other  digital  and  analog  modulation  types  must  be  investigated 
before  the  algorithms  can  be  fiilly  tested. 

5.  Other  Signal  Models 

Models  D  and  E  more  accurately  reflect  the  situations  most  often  encountered  by 
systems  attempting  to  measure  TDOA  and  FDOA  for  SOIs.  Algorithms  derived  directly 
from  these  models  could  prove  to  perform  even  better  than  those  derived  under  gimpiifi<>it 
conditions.  While  seemingly  mathematically  intractable,  an  investigation  into  this  method 
seems  worthwhile  given  the  possible  benefits. 

6.  Lower  Bound  Extensions 

The  extension  of  the  CRLB  to  the  cyclostationary  TDOA  case  has  been 
accomplished  and  reported  by  Schell  in  [6].  According  to  Gardner  [23],  this  extension 
would  not  apply  to  the  TDOA/FDOA  case  explored  here.  In  addition,  the  reportedly 
tighter  Barankin  and  Ziv-Zakai  bounds  have  yet  to  be  investigated  for  application  to  the 
cyclostationary  case.  Work  on  these  LB  extensions  could  be  highly  valuable  to  furthering 
the  understanding  of  the  beneficial  effects  of  modeling  signals  as  cyclostationary. 

7.  Other  Optunizatioii  Techniques 

Using  a  different  optimization  technique,  an  investigator  might  discover  a  more 
computationally  efficient  algorithm  that  performs  equally  or  better  thanjthe  approximated 
ML  algorithms  developed  here.  A  method  such  as  Maximum  Entropy  might  be  employed 
as  a  different  approach  with  positive  results. 
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APPENDIX 


Al.  SIGNAL  OF  INTEREST  PARAMETERS 


San^Iing  Frequency: 

/s=lMHz 

Carrier  Frequency: 

Keying  Rate: 

Bandwidth: 

5soi=(%6)-^» 

TDOA: 

Dsoi  =  50 

FDOA: 

^•‘  =  (^56)-^* 

A2.  WIDEBAND  SIGNAL  NOT  OF  INTEREST  PARAMETERS 


San^ling  Frequency: 

/s=lMHz 

Carrier  Frequency: 

Keying  Rate: 

/k  =  0.1/s 

Bandwidth: 

Bsoi  =  0.3/s 

TDOA: 

Dsoi  =  58  ^ 

FDOA: 

A3.  NARROWBAND  SIGNAL  NOT  OF  INTEREST  PARAMETERS 


Sampling  Frequency: 
Carrier  Frequency: 

Keying  Rate: 
Bandwidth: 

TDOA: 

FDOA: 


/s=lMHz 

A  =  0.025 /s 
Bsoi  =  0.075 fs 
Dsoi  =  58  ^ 
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